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Thesis Abstract 



In the first part, we investigate the effect of long range particle exchange in ideal bosonic- 
chains. We establish that by using the Heisenberg formalism along with matrix product 
state representation we can study the evolution as well as the ground state of bosonic 
arrangements while including terms beyond next-neighbour hopping. The method is 
then applied to analyse the quench dynamics of condensates in a trapping potential and 
also to study the emergence of entanglement as a result of collision in boson chains. 
In the second part, we study the ground state as well as the dynamics of ID boson- 
arrangements with local repulsive interactions and nearest-neighbour exchange using 
numerical techniques based on time evolving block decimation (TEBD). We focus on 
the development of quantum correlations between the terminal places of these arrange- 
ments. We find that long-range entanglement in the ground state arises as a result of 
intense boson tunnelling taking place across the whole chain in systems with appropriate 
hopping coefficients. Additionally, we identify the perturbations necessary to increase 
the entanglement between the end sites above their ground state values. In the final 
part, we study the wave function of a kicked condensate using a perturbative approach 
and compare the results obtained in this way with numerical simulations. 
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Frequently Used Abbreviations 



EEE = End-to-End Entanglement 

BH = Bosc-Hubbard 

PTH = Perfect Transmission Hopping 

CH = Constant Hopping 

MPS = Matrix Product States 

TEBD = Time Evolving Block Decimation 

Note About Units 

Throughout this work we measure energy in units of the recoil energy En, an energy 
reference very common in optical lattice experiments. In most cases, we explicitly indi- 
cate the energy units, otherwise, it should be assumed that energy in being measured 
in terms of the recoil energy. Similarly, we use the dimensionless parameter as a 
measure of time. We have chosen not to give units to variables that represent imaginary 
time, because such variables do not have a direct physical meaning. They are given in 
arbitrary units. 

Note About Graphs 

In all our simulations of bosonic chains, we always consider chains of unit filling, that 
is, the number of bosons is equal to the number of sites. 
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Introduction 



Weakly interacting systems can be studied using few-component models, which describe 
the physics of a small number of particles isolated from their environment. Such an 
approach has been successful in reproducing, thoroughly or partially, a large variety 
of physical phenomena studied since the establishment of quantum mechanics over the 
past century. However, particles in real systems, specially strongly correlated systems, 
interact with each other and develop a long-scale coherence that causes deviations from 
the predictions of simple models. As a result, understanding the physics of highly 
correlated systems has become the focus of contemporary physics. 

Bosons and fcrmions obey different statistical properties that determine the behaviour 
of compound systems, specifically, bosons can occupy the same quantum level while 
fcrmions cannot. Nowadays, the tremendous sophistication of cooling techniques in op- 
tical lattices allows a closed-form study of atomic and molecular systems in combination 
with optical interactions. As a result, physicists have been able to probe not only single 
particle physics in weakly interacting phases, but also the arising and taking over of 
highly correlated states of matter. This has prompted a major interest in strongly in- 
teracting systems whenever the resources to observe many-body effects under controlled 
circumstances are now available using state-of-the-art technology, which, at the same 
time, has revolutionized the way as scientists approach both theory and experiment. 
Indeed, while a couple of decades ago the characteristics of the sample under study 
depended almost entirely on its inherent physical composition, today it is possible to 
create samples with desired properties and characteristics using optical lattices. One 
of the most outstanding achievements in experimental physics occurred in 1995 with 
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the observation of boson condensation in optical latices of 87 Rb and 23 Na reported by 
Anderson et al. PQ and Davis et al. [5| respectively. Such observations have since been 
the subject of intense theoretical and experimental investigation. Certainly, the reason 
for this growing interest is twofold. On the one hand, people are interested in practi- 
cal applications, while on the other hand, there is a compelling desire to scrutinise the 
canonical framework that sustains contemporary physics. It is because of these reasons, 
and others which we will point out further ahead, that we have opted for concentrating 
on bosonic models. 

Systems of bosons display several unique characteristics, but it was the observation of 
superfluidity in 4 He that ultimately triggered the scientific desire to understand the 
physics behind the many-body effects of bosonic systems. Ever since, superfluidity has 
been explained in terms of the Bose- Hubbard (BH) model, which describes a system of 
bosons with hopping and repulsion. An alternative version of the BH model, known as 
the Hubbard model, can be used to study fermions, but here we focus on the BH model 
unless otherwise stated. Notably, it has been recently shown by Ho et al. [3] that the 
phases of the Hubbard model can be simulated by a version of the same Hamiltonian 
with attractive, rather than repulsive interactions. On the other hand, it was Fisher 
et al. [4. in 1989 who first unified the fragmented knowledge available at the time 
and established a consistent theoretical framework for the BH model. In that work, 
the BH model is studied both in the absence and presence of disorder, and parametric 
phase diagrams depicting the different phases of the model are discussed. The standard 
BH model shows two basic phases, the Mott insulator and the superfluid. The first is 
characterized, among other features, by the existence of an energy gap. The transition 
from insulator to superfluid is found to be mean field in character and universality 
properties are also discussed. Conversely, it is found that in the presence of disorder a 
new phase in between the insulator and the superfluid exists. This is called the Bosc- 
glass and is similar to the Mott insulator, but has no gap. Importantly, it is argued 
that in the presence of disorder the transition to the superfluid is always from the 
Bose-gas, and never directly from the Mott insulator. This work settled the basis for 
future approaches to the BH model, but the lack of numerical methods and additional 
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experimental applications prevented further advancements for a short time, even though 
additional investigations were carried out in subsequent years. It was not until 1999 with 
the article of Jaksch el at. [5] that an experimental proposal to verify the BH model 
in optical lattices was published. It was shown that both the Mott insulator as well 
as the superfluid regimes were reachable in optical lattice experiments. Three years 
later the transition in a gas of 87 i?6 was observed in the experiment of Greiner et al. 
[3] , using cooling techniques previously employed in Bose-Einstein condensation. In this 
experiment the phases were identified from absorption images after ballistic expansion of 
atoms. Certainly, while a gas in a Mott insulator phase projects a Gaussian distribution 
with non- visible coherent features, images from the superfluid gas display fringes, which 
arc interpreted as a signature of Bosc-condcnsation of momentum wave-functions. A 
variety of experiments has been taking place ever since the presentation of this pioneering 
work. Here we would like to mention the work by Stoferle et al. [7], where the phases 
of the gas are probed spectroscopically. Once the gas has been trapped and cooled 
in a magneto optical trap (MOT), an optical excitation is sent through the sample in 
the form of a shaking optical potential. The absorption profile is worked out from the 
absorption images after ballistic expansion. In this way, the Mott insulator can be 
identified from the peaks of the absorption profile, which indicate a coincidence between 
the energy of the optical excitation and the insulator gap. Curiously, absorption images 
from the superfluid phase indicate optical absorption stronger than the one observed in 
the Mott phase. This unexpected behaviour does not match the phase diagram of the 
Hubbard-model, since the superfluid is essentially gapless. This may indicate that strong 
many-body effects supersede particle tunnelling in sections of the phase diagram where 
the superfluid is due to exist. Likewise, the marked absorption profile has been recently 
used to probe electromagnetically induced transparency (EIT) in Mott insulators, as 
reported by Schnorrberger et al. in reference [5]. Another interesting effect observed in 
optical lattices is the condensation of fermions in the form of Cooper pairs reported by 
Regal et al.\§\ and Bourdel et al. [TU] in 40 K and 6 L respectively. In these experiments 
the interaction among particles is varied using a Feshbach resonance, which can be 
induced using a magnetic field applied directly on the sample. In the repulsive regime, 
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the fermionic atoms couple in dimers, forming weakly bounded bosonic-molecules that 
can undergo Bosc-Einstcin condensation. In the attractive regime, on the other hand, 
fcrmions couple in Cooper pairs, and then the pairs condense in the lowest energy 
level. In this latter case experimental detection is challenging as the fermionic nature of 
particles prevents ballistic expansion, therefore alternative techniques are implemented. 
In the same way, equally revealing experiments in optical lattices have been reported over 
the past years, a few of which we include in our references [TTJ [T2J H31 [HI [HI HI H3 HB] ■ 
Simultaneously to the development of cooling techniques and the increasing experimental 
efforts in optical lattices, there have been noticeable advancements regarding many-body 
numerical methods. As it is well known, the description of real quantum systems require 
exponentially large resources. As a result, the use of numerical approaches becomes es- 
sential. One of the main breakthroughs came in 1992 with the work of White [19j [20] . 
and the introduction of the density matrix renormalization group (DMRG) method to 
calculate the ground state of many-body systems. Since its introduction, DMRG has 
been extensively applied, with diverse emphases and enhancements, to the study of 
many-particle configurations in multiple scenarios and it is considered one of the most 
efficient and reliable numerical methods to date. Moreover, following the ideas under- 
lying DMRG such as the description of the system using matrix product states (MPS), 
another method was proposed by Vidal in 2003 to simulate both real and imaginary time 
evolution of slightly entangled systems [H| . The method was later introduced as time 
evolving block decimation (TEBD) [35]. These works encouraged further research in the 
area in subsequent years such as implementations in infinity systems, known as iTEBD 
[2"5] , the use of disentanglers to increase the simulation efficiency and applications 
in two dimensions using the so called multiscale entanglement renormalization ansatz 
(MERA) [35], among other contributions by Vidal's group. Equally important studies 
have been carried out by the group of Verstraete et al. ([M] and references therein), who 
have utilized MPS to describe mixed states. They also have applied MPS to simulate 
the master equation and have introduced collateral methodologies based in what they 
call matrix product operators (MPO), in contrast to MPS. Similarly, in the work by 
Hartmann et al. [271 128j it has been shown that in certain circumstances the use of 
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Heisenberg operators has advantages over the usual approach. Additionally, a method 
based in DMRG that can be used to simulate time evolution of many-body systems, 
therefore known as tDMRG, was introduced by White and Feiguin shortly after 
Vidal's seminal paper. 

This development in numerical methods along with the increasing availability of com- 
puting technology has provided the tools to explore highly correlated systems. However, 
the application of such numerical methods to relevant physical models is by no means 
straightforward. In fact, each problem possesses its own set of complications and handi- 
caps. The first works regarding the use MPS-alikc methods in highly correlated systems 
often considered, among other systems, the Hubbard model (as for example in [30]). 
which needs a supporting space smaller than the BH model. Not long ago a complete 
numerical study of the BH model in one dimension including next-neighbour interac- 
tions was undertaken by Kuhner et al. [3TI [32] using DMRG, but the first application of 
TEBD in BH chains came with the work of Daley et al. [33] , where the currents resulting 
from a density gradient in a ID bosonic arrangement are studied in the presence of an 
impurity in the centre of the chain. The impurity works as a switch (transistor) that 
can be used to control the flux of particles in a process that resembles the phase inter- 
ference effect underlying EIT. The authors implemented an enhanced version of TEBD 
that uses the conservation of the total number of bosons to improve the efficiency of the 
simulation. The currents are characterized in terms of the system parameters and very 
interesting results are shown, although the authors report a number of sensitive issues 
regarding the behaviour of TEBD, which seemed to reproduce inconsistent results under 
specific circumstances. A similar approach was explored in the paper of Hartmann and 
Plenio |34| , but this time the currents resulting from a difference in the phases of two ad- 
jacent chains are the focus of study. Namely, a chain prepared in a Mott insulator state 
is connected to a chain of equal size prepared in the supcrfluid state. Particles migrate 
from the insulator towards the supcrfluid generating a bosonic current that is simulated 
using a symmetry-enhanced TEBD. Similarly, the method was used by Mishmash et al. 
[551 131)] to simulate the evolution of dark solitons in a chain of ultracold atoms. Inter- 
estingly, simulations show that soliton waves lose coherence and eventually vanish as a 
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result of non-linear effects induced by the BH Hamiltonian. In a different investigation, 
Muth et al. [37] analyse the phase diagram of the BH model in the presence of disor- 
der using both TEBD and iTEBD. TEBD has also been utilized by Mathey et al. in 
reference [35] to identify supersolid phases in ID boson-mixtures. Additionally, TEBD 
has been used to simulate the response of a bosonic system to a sudden displacement 
of the confinement potential in the letter by Danshita and Clark J3S]. In this work a 
first principle approach to the experiment of Fertig et al. |14j was proved successful. A 
similar paper by Montangero et al. (40] employed tDMRG to explain the observations 
of the same experiment. 

As can be seen, TEBD has turned out to be a very useful tool in the study of bosonic 
systems with a growing interest from the scientific community to apply the method to 
diverse situations and scenarios. Motivated by such enthusiasm, in this work we have ap- 
plied the method and its underlying ideas to the study of entanglement in boson chains. 
As it is now well known, entanglement is the main resource of quantum information 
processing (QIP) [41] . and as such it has received much attention and analysis. In the 
context of many-body problems, entanglement has been studied extensively in arrange- 
ments of two-level systems such as spins, qubits and fermions [4"2l 14"31 l4"4l 1431 05] (for 
a complete review of entanglement in many body systems, including bosonic systems, 
see reference [H]), where the size of the local Hilbert space of each site is bounded by 
a small integer. Boson chains, however, do not necessarily offer the same advantage, 
as the associative nature of bosons demands a broader spectrum of states in order to 
fully characterize the quantum state. Studies of entanglement in bosonic systems have 
focused on diverse kinds of entanglement [35] HH] , but here we focus specially on the en- 
tanglement shared between the end sites of the chain. It has been already demonstrated 
by Campus- Venuti et al. [42] and Eisert et al. [SH] that distant places of a quantum chain 
can be entangled without the need for a direct interaction among them. In addition, 
the problem of calculating the entanglement between the terminals of a boson chain has 
been attacked before in the works of Romero-Isart et al. [51] , where the dynamics is 
reducible to a single-particle propagation, and Plenio et al. [52] using Gaussian states. 
Conversely, the situations analysed here include a number of effects that do not allow 
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the reduction of the problem as mentioned just before, and therefore the use of TEBD 
becomes essential. We found that in the standard BH model the entanglement between 
distant places of the chain decreases as the chain-length augments, but we also found 
that the same entanglement can be made to increase using a slightly modified version of 
the Hamiltonian with variable, rather than fixed hopping constants. Our results are ex- 
plained in terms of the tunnelling profile displayed by the particles along the chain. We 
argue that in the case where a form of strong and resilient entanglement arises, the chain 
undergoes a special kind of fluidity enhancement in which particle tunnelling takes place 
across the whole length of the chain and not only inside localized clusters of sites. We 
hope that our results provide a basis upon which further advancements could be made 
in the same direction. In this document we also discuss how to implement alternative 
numerical methods based on MPS that can be used in a variety of situations where 
the Hamiltonian is sufficiently regular to allow an explicit solution of the Heisenbcrg 
equations of motion for the operators. The methodology proposed is then employed to 
simulate the propagation of bosons with emphasis on the amount of entanglement gen- 
erated during the evolution. As a complement to our studies of the Hubbard model we 
also present an analytical development of a system consisting of atoms under the action 
of a periodic kicking. In this part we use the Gross-Pitaevskii equation to generate the 
dynamics of atoms, which are considered as a single cloud and not individually as in 
the BH model. In a sense, we can say that this study has been inspired by the works of 
Zhang et al. [53j and Liu et al. [54j . where basically the same phenomenon is analysed 
following the experimental realization of Moore et al. [SH] in ultracold atoms. 

From our point of view, finding the conditions for the emergence of entanglement be- 
tween distant places of a boson chain is an important contribution, not only because 
long range correlations are crucial in quantum information protocols, but also because 
such correlations give insight about the system phenomenology. Another merit of the 
work is the difficulty associated with some of the numerics that we present below, espe- 
cially those regarding entanglement. This is because long range entanglement, especially 
entanglement between extreme places, inevitably involves all the intermediate degrees 
of freedom in between the boundaries, and so one needs to synchronise a large number 
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of correlated processes that derive from our use of MPS instead of a number-of-particles 
basis. On the other hand, we consider that the method that we introduce further ahead 
will prove useful in the study of dynamical models as it does not depend so heavily on 
the amount of entanglement in the system as the same TEBD or DMRG, which makes 
it suitable to perform long-time simulations, although our method does not apply to the 
same wide spectrum of problems as the previously mentioned do. 

This document is organized as follows, in chapter Q] we include a review of some of the 
concepts and methods that are to be used in subsequent chapters such as entanglement 
and matrix product states. Chapter [2] sketches our numerical approach and discusses 
several programming issues that proved sensitive while coding our algorithms. Chapter 
[3] deals with the application of MPS in circumstances where the regularities of the 
Hamiltonian allow an elegant application of this representation. Chapter [4] introduces 
the BH Hamiltonian and shows how TEBD works specifically for this model. In chapter 
Owe apply the method to find the ground state as well as the dynamics of the BH model 
under different circumstances. A number of schemes are proposed and analysed. We 
then proceed to present an analysis of the dynamics of cold atoms driven by a periodic 
kicking in chapter [51 After this, we summarize and present our conclusions. 
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Chapter 1 



Preliminary concepts 

1.1 Schmidt decomposition 

Given a pure quantum state of a system made up of many individual components, it is 
possible to write the quantum ket of the system as a sum of product of states in the 
following manner [41 j . 

IV) = ( L1 ) 

i 

where kcts \i>i) and form orthonormal sets of vectors corresponding to complemen- 
tary subspaces. Namely, 



(vi\vj) = 6\, 

{ Vi \ N ) = 0, (1.2) 

where 6f is the Kronecker delta. These kets are known as the Schmidt vectors. Similarly, 
the set of Aj are the Schmidt coefficients. The Schmidt coefficients are positive real 
numbers that satisfy the condition, 
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J> a = l. (1.3) 

i 

The process of writing the state as in equation is known as the Schmidt decom- 

position. Very important consequences follow from the Schmidt decomposition. For 
instance, it can be shown that the reduced density matrices describing the states of 
complementary subsystems share the same set of eigenvalues, which are equal to the 
square of the Schmidt coefficients. Similarly, there are two important characteristics 
that we want to remark. First, the Schmidt decomposition depends on how the origi- 
nal system is divided into complementary subspaces. This division may or may not be 
related to the actual geometry of the system. Second, the Schmidt decomposition for a 
given partition is in general not unique. For instance, given a particular set of Schmidt 
vectors we can obtain a different set of Schmidt vectors for the same partition by ap- 
plying unitary operations on the original vectors. This however does not produce any 
change on the coefficients, which remain as positive numbers. In a sense, the topology 
of the decomposition is preserved after unitary operations, but this is only true when 
such operations take place in the subspaces that support the Schmidt vectors. One way 
of getting the Schmidt vectors is by simple inspection. Obviously, this is not at all 
practical when we are dealing with intricate quantum states. The standard method to 
get Schmidt vectors is by diagonalizing the reduced density matrices corresponding to 
the subspaces defined by the partition. As we will see, this property is at the heart of 
the numerical techniques that we will use to study our models. 



1.2 Entanglement characterization 

Even though the fundamental concepts of quantum mechanics were already well estab- 
lished and accepted by the scientific community by 1930, it has not been until recent 
times that the concept of entanglement has started to receive attention. Among other 
reasons, it is because nowadays people arc quite interested in what aspects of the physical 
systems are purely "quantum" . Entanglement is a characteristic associated exclusively 
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with quantum states. Classical representations of physical systems do not contain any 
form of entanglement whatsoever. It is for this reason that entanglement provides a 
measure of how efficient a physical system can be and how much of the state provides 
quantum resources that can be potentially used. This is why quantum entanglement is 
now a new branch of physics with a promising future. Quantum entanglement is the 
main resource of several highly efficient tasks proposed in QIP such as supcrdense cod- 
ing, quantum state telcportation and quantum cryptography. Quantum entanglement is 
also at the heart of the quantum computer. In many senses, the quantum computer can 
be more efficient than its classical parallel. Let us take for example one of the simplest 
tasks of a computer: generate random numbers. Notably, this apparently simple oper- 
ation carries some difficulties for a classical computer, as such a machine is essentially 
deterministic. Usually, random distributions are generated using recursive algorithms 
that always introduce deviations. For a quantum computer, on the other hand, the 
generation of a random distribution would result naturally by performing measurements 
over an equivalent state superposition. This simple example illustrates the usefulness of 
quantum states in a practical scenario. Furthermore, it is the concept of entanglement 
what captures the degree of utility of a quantum state. As a matter of fact, entangle- 
ment has proved to be a rich field of theoretical investigation from which very interesting 
results have been derived [561 1571 158] . The progress in experimental physics has been 
interesting but not equally dynamic, although a number of experiments involving en- 
tangled states have been carried out with relative success [T7]. So far, entanglement 
appears to be consistent with experimental observations, but practical implementations 
using entanglement as a resource remain challenging |18j . 

Entanglement in a pure bipartite system can be consistently characterized from the 
reduced density matrix of any of the component subsystems. Let us call pa the reduced 
density matrix of subsystem A, ps with an analogous meaning, then entanglement 
between subsystems A and B is given by the von Neumann entropy, 

S = -tr(p A log 2 p A ) = -tr(p B log 2 p B )- (1.4) 
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Furthermore, S can be calculated directly from the eigenvalues of either density matrix, 



It can be verified that S = for a separable state. In fact, when the state is separable 
the reduced density matrix contains only a single eigenvalue equal to one so that the 
logarithmic function in (|1.5[) causes the whole expression to vanish. It can also be 
verified that the maximum value displayed by S is log 2 [d], where d is the dimension 
of the smallest subsystem. Von Neumann entropy provides a reliable estimation of the 
amount of entanglement shared between subsystems A and B as long as the whole AB 
system remains in a pure state. In most cases von Neumann entropy is an operational 
criterion, which means it can be calculated from the expression that gives the quantum 
state. Such is the case when the state is given in terms of a discrete basis. Then, in order 
to get the entanglement, the reduced density matrix must be found and S is computed 
from the eigenvalues of such matrix. This procedure can be considerably difficult to 
apply when the state is given in a continuous basis. In this case the equivalent of finding 
the eigenvalues of the reduced matrix corresponds to solving a second order differential 
equation. Therefore, it is more convenient to use a criterion such as the trace of p A 
which can be obtained from direct integration. The amount of entanglement is in this 
way given by how much the trace deviates from one. Similarly, more entanglement 
criteria can be worked out, but S is well established as the most consistent measure of 
entanglement for bipartite pure states. Among the conditions that a good entanglement 
measure, say E, should satisfy, we find [58l 159] . 

• E must be positive 

• E must be zero for separable states 

• E must be invariant under unitary transformations performed on subsystems A 
and B. Hence, we can understand why S depends entirely on the eigenvalues of 
the reduced matrices, precisely because the eigenvalues are invariant under unitary 
transformations. Additionally, because reduced density matrices of subsystems 




(1.5) 
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that correspond to complementary partitions of a pure general state share the 
same set of eigenvalues, S is independent of which reduced matrix one chooses to 
make the calculation. This is not the case, however, when the state of the whole 
system is mixed, since the eigenvalues of reduced matrices obtained from reducing 
a bigger matrix may be different. 

• E must not increase under local unitary operations and classical communications 
(LOCC). Actually, it can be shown that this property implies the previous one 
in the absence of classical communications. Such classical communications make 
reference to the sharing of information between parts A and B using classical 
resources such as telephones or computer networks. 

• There must be maximally entangled states. This implies that E establishes an 
ordering of elements in the Hilbert space. This condition cannot be fully incorpo- 
rated for multi-particle entanglement-measures and therefore the statement should 
be considered only in the bipartite realm. 

The fact that reduced density matrices of mixed states do not share the same set of eigen- 
values as in the case of pure states prevents a direct generalization of results such as 
equation (|1.4[) . One of the main breakthroughs in the subject of entanglement character- 
ization is due to Peres and members of the Horodccki family [601 161] who simultaneously 
came up with the idea of partial transpose, which we now present. A general separable 
mixed state of systems A and B can be written as, 



So, if we swap matrices pf by their corresponding transpose matrices, the expression 
above would still be a valid density matrix since the transpose matrices of pf are density 
matrices themselves. As a consequence, the new density matrix describing the whole 
system is a positive operator, that is, all its eigenvalues are positive. This operation 
can be generalized to the case when matrix p is not separable. In such circumstance 
transposing subsystem B would involve swapping the indices of the big density matrix, 




(1.6) 
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namely, 



(1.7) 



where we have used 77 as the partial transpose of matrix p. Crucially, in this case we 
cannot argue that r) is a valid density matrix and a positive operator, therefore, showing 
that 77 is not positive provides evidence that the state is entangled. Consequently, in 
order to know whether the state is entangled, it suffices to get the partial transpose and 
see if any of its eigenvalues turn out to be negative. On the other hand, it is also good to 
comment that this criterion does constitute a necessary rather than sufficient condition 
for entanglement. It could be that A and B are entangled albeit rj having a positive 
spectrum. It is possible to take one step forward and formulate a quantitative expression 
for entanglement from the partial transpose criterion. Vidal and Werner |62| propose 
Log-negativity, which bounds the amount of distillable entanglement in p. By definition, 
distillable entanglement makes reference to the pure state entanglement that can be 
extracted from p. Log-negativity is an extension of the simpler negativity, the sum of 
the negative eigenvalues of the transpose. Formally, Log-negativity can be written as, 



where are the eigenvalues of 77. According to |62j . Log- negativity is an additive 
quantity, which makes it a very useful measure. Also, the fact that it is given in terms 
of the logarithm function allows a comparison with von Neumann entropy in systems 
of zero mixture. In this case, it has been shown that Log-negativity provides a greater 
estimation of entanglement than von Neumann entropy. We want to emphasise that 
von Neumann entropy can only be applied to measure the entanglement of pure states. 
In this thesis, the entanglement of mixed states will be calculated using the definition 
Log-negativity presented above. In a sense, Log-negativity is a quantification of the 
Peres-Horodecki criterion. One positive consequence of this relation, is that when a 
quantum state is shown to be entangled according to the Log-negativity criterion, then 
it can be shown that the entanglement contained in such state can be distilled, that is, it 




(1.8) 
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can be transformed into pure state entanglement through a distillation process. As it has 
been mentioned before, Log-negativity sets an upper bound for the amount of distillable 
entanglement of a given state. Hence, we can say this entanglement measure establishes 
the degree of usefulness of a quantum state. Sometimes it is interesting to see the 
relation between entanglement and correlations. To this end, it is important to highlight 
that entanglement is a property associated exclusively with the state, while correlations 
require the intervention of an observable. It is widely accepted that entanglement is a 
resource more exotic than correlations. There can be correlations without entanglement 
but there cannot be entanglement without correlations. In the models studied in this 
thesis, correlations play an important role in defining the phases of the system. We 
explore the development of entanglement in situations where correlations are known to 
induce highly intricate states. As we will see, this can derive in very entangled states 
that we try to characterize from the physical features displayed by the system. 



1.3 Matrix product state representation 

Initially, when quantum mechanics began to be used as an accurate theory capable of 
delivering insight in atomic systems, fixed bases sufficed to provide an operational plat- 
form to perform, most of the time, analytical calculations. Even when Dirac introduced 
the interaction picture, a method in which basis kets and operators evolve according to 
the dynamics generated by the integrable part of the Hamiltonian while the quantum 
state evolves according to a non-diagonal term, it was intended to be used mostly as a 
complement of the existing quantum pictures, namely, Hciscnberg's and Schrodingcr's. 
The idea of a dynamical basis has been brought recently, partially as a result of the in- 
sight obtained in the process of understanding entanglement. Let us think of dynamics 
as a unitary operation that is applied on the initial state. If we want to make things easy 
and use a basis in which the evolved state remains simple at all times we would rather 
employ a basis that does not change, or changes little at least, when unitary operations 
are applied on the quantum state. For someone who is familiar with the formalism 
of entanglement, the idea of Schmidt vectors is likely to come to mind in this specific 
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In> |i> Iv> 




Figure 1.1: Example of MPS. and \v) are Schmidt vectors from the decompositions 
at the sides of place i, 

situation. Indeed, Schmidt vectors do not change when unitary operations are applied 
on the system. As a consequence, entanglement between complementary subsystems is 
invariant under local unitary operations. Nevertheless, the evolution operator acts glob- 
ally and therefore induce changes on every Schmidt vector along the system. However, 
when the evolution operator can be split, at least approximately, into non-overlapping 
semi-local operators, the idea of using Schmidt vectors as basis vectors becomes feasible. 
Such is the case for those spin and boson chains where hopping takes place only among 
next neighbours. In such a scenario every Schmidt decomposition determines a splitting 
of the chain. Schmidt vectors describe the state of a subset of the chain. In general, the 
Hamiltonian can be written as, 

H = J2 (Kk+i + B k ) , (1.9) 

fc 

while the evolution operator is given by 

U = e~ mA « Y[e- m ( Ak - k+1+ii «) =Y[U k (1.10) 

k fc 

In this way, unitary operations applied on pairs of neighbour sites do not affect all the 
Schmidt decompositions. This suggests that Schmidt vectors can form a basis suitable 
for state description, one in which the process of state updating is efficient. To see how 
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Schmidt vectors can be used, suppose that the chain is in a pure state. To see how the 
state can be written in matrix product state (MPS) representation, let us focus on one 
site of the chain. Initially, in order to describe the state, we use, on the one hand, a 
local representation for such a site, say \i), given in a standard basis, for instance spin 
orientation or occupation number, on the other hand, to describe the state outside this 
site, we use Schmidt vectors obtained from splitting the system to both sides of the 
place in consideration (figure [TTTJ) . So, making explicit reference to a local basis at site 
n the quantum state reads, 



M =E^'- 1 lAWr| I f| M )i 1 -- 1 ]|i)[''V)K^]. (1.11) 

///// 

In this way written, the quantum state appears as a superposition of states in the basis 
of the Schmidt vectors and the local basis of site n. We have deliberately specified 
the Hilbert spaces of every ket as a superscript. Tensor F^u contains the components 
that describe the state in the new basis while coefficients X^ v are just the Schmidt 
coefficients, necessary in this case to guarantee that vectors and \v) are normalized. 

One important property of this tensor decomposition is that it is actually possible to 
establish a relation among the T's and the Schmidt vectors via induction. Indeed, vectors 
\v) can be expanded in terms of a basis and Schmidt vectors \^)i n + 2 ' N i in this 

way, 



^[n+hN] = ^ r |«+lb A [n+l]|^[n+l]|^[n+2,W] 

In fact, according to [3T], and \v) can be written in terms of the state tensors 
associated with the subsets [1, ...,n — 1] and [n + 1, N] and the standard basis in the 
following form, 

l^ 1 '- 11 = E f E r^uW...r["-^- 1 K l! ... ) i n - 1 )], (Lis) 



and, 
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W)l n + i, N] = Y, r^ 1 ^...^"^^!^,....^ , (i.i4) 

7,...,w \ i n+1 ,...,ijv / 

where N is the number of sites in the chain. Additional relations can be derived, in 
particular we would like to mention that the standard coefficients can be put in terms 
of the coefficients of the decomposition in the following way, 

c U,S2,...,iJv — ^ 1 la A a 1 a/3 1 "1 1 l 1 - 10 ^ 

a,/3,...,u; 

These coefficients enable us to write the quantum state back in the conventional basis, 

\lp) = c ii,i2,...,twl i l)*2, -,ijv), (1-16) 

n,«2, ■■■,'» 

so that the representation of the state in terms of Schmidt coefficients and tensors, to 
which we refer to as canonical decomposition, fully characterizes the quantum state. 
The canonical decomposition is very convenient to study numerically the time evolution 
of the chain. Every time that a semi-local unitary operation as in equation (|1.10[) is 
applied, the canonical representation is to be updated only for the elements involved 
directly with the transformation. For example, if a unitary transformation is applied 
on sites 4 and 5, we can see that for instance r[Jj, 47 does not need to be updated, since 
the topology of the decompositions associated with the tensor, namely those in which 
|^[!,6] an( } |^[8,jv] are involved , is not affected by the operation. For the case of vectors 
l/x)! 1 ' 6 ] particularly, a unitary operation on sites 4 and 5 induce a change on every vector, 
but the new vectors are Schmidt vectors of the evolved state with the same Schmidt 
coefficients. So, the new Schmidt decomposition preserves exactly the same configuration 
of elements with exactly the same coefficients of the initial decomposition. Therefore, no 
change takes place on the elements of the decomposition in places far away from where 
the semi-local unitary transformation operates. Following the same analysis one can 
establish that a unitary operation on two neighbour sites n and n + 1 generates changes 
only in r^, rl n+1 l and A^. This property constitutes the basis of efficient simulation: 
as the cost of updating the state depends only on the local characteristics of the system, 
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numerical routines of low memory consumption and fast execution can be implemented 
using the canonical decomposition. Crucially, the factor that determines the speed of the 
simulation is the number of Schmidt vectors in the single value decomposition. Systems 
with few vectors can be updated through relatively few computational steps 
Moreover, the canonical decomposition is also known as matrix product states (MPS) 
representation, and extensive documentation can be found under this denomination |26j . 
To see how the canonical decomposition can be updated after a semi-local unitary matrix 
operates, let us first focus on the simplest case where a one-site unitary operation, W n \ 
is applied on an arbitrary site n. In this case we can use the expression in equation 
(jl.lip to apply the transformation directly, 

i^> = t)-wi^> = x)S) A ^" llA ^ ] ^' ,[ " lr ^*iA*> I1 '"~ 11 i*> w i' / > [ " +1 ' JV1 » ( L17 ) 

fJ,u i 

so that the effect of the unitary matrix is reflected only in the coefficients T<- n '. The 
updated canonical coefficients are identical to the originals, with only one exception, 

fe 

Furthermore, this transformation does not increase the number of elements necessary to 
describe the state and the updating procedure can be coded easily. 
Let us now assume that a two-site transformation, is applied. In order to carry 

out such operation, the state must be written with explicit reference to the local basis 
of sites n and n + 1. This is done by inserting (|1.12[) in which results in, 

IV') =EE A A r ^ <A *' r 5 +11J V>|i>li>IO- (1-19) 
Therefore, the updated state reads, 

i/) = i>["<" +1 ^) - EE A ^ A ^?' r ^ A - r S +11 ^)h'}ii>ie>- (i.2o) 

ftv£ ij 

As discussed before, this transformation induces changes in tensors r|Tj\ rJjV 1 " 1 ^ and 
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\ u as a result of the changes generated in the Schmidt decompositions of the splitting 
[l,n—l] : [ti, iV] and [l,n] : [n + l,N]. Consequently, the next step consists in finding 
the reduced density matrix of subsystem [n + 1, N], from which new updated Schmidt 
vectors can be obtained as eigenvectors: 



As a result of the reduction, wc are left with a density matrix spanned by the local 
basis of site n+1 and Schmidt vectors Crucially, the cost involved in manipulating 
such density matrix, that is, storage and diagonalization, is proportional to the number 
of Schmidt vectors in the state. In this way, the astronomical memory requirements 
associated with the standard basis can be avoided, as long as the number of Schmidt 
vectors remains small. New coefficients \' v can be identified as the square roots of the 
eigenvalues of the reduced density matrix, whose eigenvectors can be used to get the new 
tensor r'!^ 1 ^ directly from equation (|1.12j) . Finally, the new tensor T'jfJ 1 comes from 
projecting the Schmidt vectors over the whole state given by equation (|1.20p [2"Tl 125] , 

In order to exemplify how the canonical decomposition can be used to represent the 
state, let us focus on a system of three qubits. Suppose that the state of the system is 
given by, 



Now, if wc take the first qubit and consider the other two as the rest of the system, we 
can see there is only one Schmidt vector to the right, namely, \v) = |10). Then, using 
the convention introduced in equation (jl.lll) the state can be written as, 



(1.21) 



= ]H0>. 



(1.22) 



|*> = AWrW 



\Q)\v)+X [ Mi\m, 



(1.23) 



with 



1 (the Schmidt coefficient) and, 
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r [i]o = o 
1 1,1 — u ' 

r [i]i _ l 
1 1,1 — L - 

Note that the sub indices of rW make reference to the Schmidt vectors to the left and 
right of the first site. In this case there is only one vector to the right while for the 
left we imagine there is an ancillary state. Similarly, for the second site we can see the 
Schmidt vectors to the right and left are = |1) and \v) = |0) respectively. The state 
can therefore be written as, 

\d>) = AWrgV)|o>k) + AWrP\V>|i)W, (i-24) 

with, 



p[2]0 _ „ 

1 1,1 — u ' 

p[2]l _ i 
1 1,1 — L i 

and =1. Finally, the coefficients of the third site can be extrapolated from, 

|0}=A?MlV}|O)+A[M 3 iV}|l}, (1-25) 

with, 



r [3]0 _ 1 
1 1,1 — x ' 

p[3]l _ Q 
1 1,1 — u > 

and as the reader may have guessed, = In the same manner, we can obtain the 
canonical decomposition of the following less trivial state, 
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looi> + |on> _ r jo>+ji> 



v/2 I V2 

This state can be written with explicit reference to the coordinates of the first qubit 
exactly as in equation (|1.23p . but with, 



r [i]o _ i 
1 1,1 — 



1 1,1 — u > 



and, 



Similarly, the state also adopts the form of equation (jl.24|) . but with the following 
important changes, 



p [2]0 _ J_ 

M " v/2' 

r [2]l _ J_ 

and |/x) = |0) and = |1). Likewise, it can be seen that for the third site the coefficients 
of equation p. 251) are given by, 



r [2]0 _ „ 
1 1,1 — u ' 

F [2]l _ i 
1 1,1 — L ' 



and the Schmidt vector is, 
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More complex representations can be derived when there are partitions with more than 
one Schmidt vector in the decomposition. We hope that these simple examples allow the 
reader to grasp an idea about the mechanics associated with writing a quantum state 
using MPS. 



1.4 Perfect transmission hopping 

In a chain of bosons with constant chemical potential and zero repulsion the Hamiltonian 
can be written as, 



N-l 

H =22 J k{a\ +1 a k + a\a k+ i), (1-29) 
fc=l 

where a\ and ctk are the usual bosonic operators with standard commuting rules, namely, 



[fife, ai] = 0, 



= 0, 



a u a\ 



SI (1.30) 



and N is the number of sites in the chain. One question of interest in several branches 
of physics is whether it is possible to dynamically transfer an arbitrary quantum state 
from one end of the chain to the other with maximum fidelity. In other words, whether 
it is possible to have a perfect transmission channel. This problem has been extensively 
studied and here we limit ourselves to outline the main results of more specialized investi- 
gations presented elsewhere [64, 65, 66, 67, 68:. Perfect transmission was first studied in 
spin chains due to the interest prompted by the novel scheme shown in reference [43] . In 
this reference, spin chains were proposed as alternative channels to transfer quantum in- 
formation encoded in the state of the spins. In this context, it became important to know 
under which specific circumstances a spin chain could transmit a state without the state 
being corrupted by the underlying dynamics. It turned out that chains with constant 
coefficients could be used as efficient transmission channels only in chains of maximum 
3 sites. However, it was also found that chains with variable hopping coefficients could 
be made efficient if the right hopping coefficients were chosen. Such coefficients could 
be identified by noticing the parallel between the dynamics of the chain and the physics 
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of the angular momentum. In this analogy, states on the ends of the chain correspond 
to cigenstates of angular momentum with large eigenvalues, while states in the centre 
turn out to be analogous to eigenvectors associated with the smallest eigenvalues. It 
was shown that in a system governed by Hamiltonian (|1.29|) perfect transmission can be 
accomplished by choosing what we call perfect transmission hopping (PTH), 



J[™ = ^VHN~k), (1.31) 

where A is a constant that fixes the time scale (for simplicity we chose A = 2 for numerical 
simulations). This stands in contrast to the widely used constant hopping (CH), given 
simply by, 



Jk = 1- (1-32) 

When inserted in Hamiltonian (|1.29|) . PTH induces a mirror-reflection of the initial state 
with respect to the chain centre at a time, 



This property constitutes the basis of perfect state transmission between the chain ter- 
minals. Independently of how many bosons initially occupy either chain end, they all 
turn up at the opposite end in a PTH chain with all repulsion constants set to zero [52] . 
Although the dynamics of a PTH chain with no repulsion has been already studied, there 
are several research extensions of interest. On the one hand, it is important to know if 
chains with variable hopping coefficients and repulsion can display perfect transmission. 
If they cannot, it would be interesting to know how the repulsion hampers the state 
transmission. Similarly, if transmission with repulsion is possible, it would be impor- 
tant to know the specific circumstances under which this phenomenon can be observed. 
Works on this direction have produced very interesting results |67| . In addition to these 
research areas, here we also focus on one alternative approach. In fact, in sections ahead 
we intend to see the problem from a rather different perspective. Certainly, it is quite 
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valid to ask how PTH affects physical phenomena in which transport does not partic- 
ipate straightforwardly. From our viewpoint, this constitutes an important aspect to 
study in chains with efficient transmission, as the regularities of the Hamiltonian that 
lead to perfect transmission could give rise to interesting coherent processes. Similarly, 
it should be mentioned that very promising experimental proposals to realize PTH in 
optical lattices have been discussed in reference [68) by the same group that proposed 
the realization of the Hubbard model in optical latices. Therefore, such kind of hop- 
ping profile should be considered as a feasible alternative and not just as a theoretical 
idealization. 
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Chapter 2 



Implementation and numerics 



The first aspect to be handled when addressing the construction of a program using the 
formalism presented in previous sections, is memory allocation. It is indeed possible 
to implement the algorithm using standard programming tools, that is, vectors and 
matrices with permanent memory attributes. However, the way in which the elements 
of the canonical decomposition depend on the Schmidt vectors makes the canonical 
tensors vary in size according to the number of Schmidt vectors in the decomposition. 
Therefore, it is preferable to adopt a programming style in which the dynamical nature of 
the algorithm can be handled more appropriately FORTRAN 95 offers several tools that 
can be conveniently adjusted to approach the dynamical nature of the method. On the 
one hand, it allows one to actually define the geometry of the objects employed to store 
data. This means that in addition to vectors and matrices, tensors of whatever shape and 
dimension can be used. For instance, we can define a vector in which every component 
is a matrix, or a matrix in which every element is made of a complex number and a 
real vector. The other tool is dynamical allocation. The elements that we use to store 
numbers are not necessarily fixed in size, instead, we can make them bigger or smaller 
according to the simulation requirements. In FORTRAN 95, dynamical allocation can 
be implemented using either allocatable objects or pointers. Pointers work very much 
as allocatable objects, but they have two useful additional properties. First, they can be 
alternatively used as aliases of other objects, and second, they can be passed from one 
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Figure 2.1: Sketch of memory allocation of LAM(i)%ld 

routine to another as arguments. Modules are another useful tool. Variables defined 
in a module can be used in any routine by just including the module in the routine's 
headlines. In our program for example, the module bskt contains the tensors of the 
canonical decomposition. In order to define an object such as A a , we must first generate 
the type that supports the object. This is done through the lines, 

type dl 

real (dp), dimension( : ) , pointer :: Id 
end type dl 

Tensor A Q , is then written as an element defined by the rules of type dl, 
type(dl), dimension(N) :: LAM 

This means that there is one pointer associated to every site of the chain. Each pointer 
is independent and can be given whatever memory one needs to allocate (figure 12. ip . 
So, if for example the Schmidt decomposition in site 5 has 3 coefficients, then in order 
to allocate that specific amount of space we write, 

allocate (LAM (5) °/.ld (3)) 

and then proceed to specify every component of LAM(5)%ld. Similarly, tensor rj™^ is 
an object that can be allocated in two dimensions, namely those associated with a and 
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/?. Every component of this allocatable arrangement, on the other hand, is a complex 
vector of fixed size. The components of such vectors correspond to a local basis, in this 
case given in terms of the label i, whose maximum value is set to be M. Therefore, we 
must utilize two types to define the variable, namely, 

type fg 

complex(dpc) , dimension(M) :: gf 
end type fg 

type dg 

type(fg), dimension( : , : ) , pointer :: gd 
end type dg 

the variable itself is defined the same way as LAM was defined before, 
type(dg), dimension(N) :: GAM 

that is, one storage unit per site. For this variable, allocation takes place only in the 
pointer section. For instance, 

allocate (GAM (2) 7.gd (3, 5)) 

and then, if we want to give numerical values to the tensor we would write, for example, 
GAM(2)7„gd(l ,3)7„gf (2) = some complex number 

In this way, no memory space is wasted through undefined memory slots. When us- 
ing memory allocation, special care must be taken in processing allocatable units. If 
a variable is allocated, then it must be deallocated if for some reason we want to al- 
locate it again, for example as the decomposition is updated, otherwise the memory 
associated with the variable will not be accessible any more and one can easily run out 
of virtual memory in a long simulation. The most challenging part of a time evolving 
block decimation (TEBD) program, is the updating routine in which a two-site unitary 
transformation is applied and some tensors must be worked out (figure (|2.2[1 ). 
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In a program that makes use of symmetries, we must first organize the vectors of the 
decomposition in groups of elements sharing the same symmetry, which at the same 
time involves finding the symmetry of the vectors themselves. Once this is done, we 
apply the updating procedure to sets of vectors with complementary symmetry, in such 
a way that one global symmetry is preserved. Finally, the symmetry associated with 
the new tensors must be found and attached in order to be able to access this important 
information in the next run of computations. The reason why we implement symmetry 
enhancement is because the cost of updating the decomposition in small subspaces is 
lower than the one of updating all the elements at once by diagonalizing a single large 
matrix. 

Similarly, in a problem where the quantum state is symmetric with respect to the centre 
of the chain, we can simulate the dynamics by updating the coefficients only on one half 
of the chain and extrapolating some coefficients in the middle. This involves using the 
fact that in a symmetric chain with an even number of sites the reduced density matrix 
of the two central sites has a degenerate spectrum and eigenvectors corresponding to 
the same eigenvalue are in fact complementary. 

One additional feature of FORTRAN 95 which hugely facilitates memory management 
are linked lists. They basically consist in variables that can be given extra allocatable 
space without destroying the information already contained in the variable. Note that 
this is not the case for pointers, since every time that a pointer is reallocated it must 
be previously deallocated and so any associated content is automatically lost. A linked 
list can be made of any valid memory unit such as vectors, matrices, pointers and even 
types. In the updating routine, for instance, we use linked lists to store the Schmidt 
vectors as we do not know in advance how many vectors conform the decomposition of 
the updated state. 

We also would like to comment about one element of the research with challenging 
numerically aspects. As it has been discussed, the calculation of entanglement in mixed 
states involves dealing with the density matrix. Here, such matrix comes from reducing 
the pure state of the whole system. In order to find the entanglement between the 
ends, or end-to-end entanglement (EEE), the reduced density matrix in the standard 
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basis must be computed using the canonical coefficients. To see how the reduction takes 
place, let us write the quantum state of the chain in the following way [21j . 

W = E |la 1 }t 1 UL 1 j|a 1 a 2 )t 2 UL1...A^-i|aiv-il) [JV1 ) (2.1) 

Oil >Ot-2,---,CiN—l 

where, 

|aa') W =E r Ll'N>- (2-2) 

i 

Consequently, the reduced density matrix looks like, 
Pin = 

E I l«i> 111 I W (aioa |a>' 2 ) ^ Ag A$ a (a 2 a 3 |a' 2 a' 3 ) [3] ... 

/ a 1 ,a 2 «JV — 1' \ 

,a' 2 ,.. ■,<*' N-l' 

...(a^ajv-ilaV^JV-i)^" 1 ^^ (2-3) 
E |lai>W(lc^i|WM {aiW _ l}Kia ^_ l} |a^_ 1 l)[ Jv l(a / Ar _il|[ w l. 

«l,a'l,a»^i,o'jv-i 

Most of the programming work is devoted to write matrix Aff aia _ N _ 1 j/ a / ia ? Jv ._ 1 \ in a way 
that can be efficiently stored and manipulated. The storage matter is related to conser- 
vation properties. In the case of Pin, conservation splits the operator in representations 
holding S ai + S aN _ 1 = Sa'-i + S a > N _^, where S a is the symmetry associated with kct 
\a). The same splitting applies for M^ aiaN1 ^ a > ia > N1 y In order to compute this ma- 
trix, one first calculates the products Aq] A^j (aio^c/ia^)' 2 ' and creates a temporary 
support matrix M/ aia2 \/ a ' lC1 ,/ 2 \. This computation can be accelerated by exploiting the 
fact that in a state like (|2.2[) a symmetry restriction bounds the indices through 

S a + S a ' + Si = S'ciobal Simmetry, (2-4) 

therefore, to compute an inner product such as, 
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(2.5) 



one just takes the indices in the kets and finds the corresponding local coordinate i 
using equation (|2.4[) . Only if the coordinates coincide the associated contribution is 
stored using linked lists. In the next run of computations, for every pair of indices 
Q!3j 3 one sweeps over the range of values of the indices 02, a' 2 and finds the product 
Aq] (0203 lava's)' 3 '- Again, only if the inner product is not zero the contribution is 
added up. After this is completed, we are left with a support matrix M{ aia3 }{ a i ia r 3 }. 
The process goes on until we finally get the matrix with the end indices. 

Improved performance can be achieved if spatial symmetry is taken into consideration. 
This is done in a very similar way than with the canonical coefficients, that is, when 
the computations come to the middle of the chain, we group up complementary repre- 
sentations from each side of the chain and perform our operations in subspaces, always 
keeping the total number of bosons fixed. 

Once the reduced density matrix for the chain terminals has been worked out, we use 
it to calculate the logarithmic negativity. This involves performing a series of rear- 
rangement operations, as for example, the explicit calculation of the partial transpose 
matrix. This process alone carries a technical issue of particular trickincss that we want 
to comment about. As we already pointed out, the conservation of the total number 
of particles allows us to split large matrices into several smaller representations that 
we can manage more easily. As a consequence, what we call reduced density matrix is 
actually a set of matrices, each one corresponding to a different quantum number. This 
can be understood by thinking that every matrix matches a state of the chain with a 
complementary quantum number. Therefore, there are as many matrices as number of 
bosons plus one, as the no-boson possibility must be accounted for too. The splitting 
is not only practical but often necessary, because the amount of information that we 
can handle using a computer is limited. Nevertheless, this splitting may be broken if 
the state is subject to non-unitary transformations. This is precisely what happens 
when we try to compute the partial transpose. As a matter of fact, the transposition 
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operation mixes subspaces with different quantum numbers. The good news is that 
the transpose matrix keeps some kind of regularity that allows us to split the matrix 
into non-interacting subspaces that we can address individually. Hence, before starting 
with the operations that determine the new matrix, we first establish the subspaces that 
conform the new splitting, and then proceed to fill every conforming matrix with the 
corresponding elements. Once this is completed, finding the eigenvalues and computing 
the logarithmic negativity from equation (jl.8p is straightforward. 

These are the numerical issues that in our opinion require special care and analysis 
before they can be efficiently implemented in a program. Even though these aspects of 
our investigation do not shed by themselves any physical insight, we have decided to 
include this discussion as a way of presenting a complete view of the problem as well as 
its collateral issues. 
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Chapter 3 

Matrix product states in ideal 
boson chains 



In a boson chain of N sites and M bosons in which particle exchange can potentially 
take place among any two places, the Hamiltonian of the system reads, 



In such a way that Rk,i represents the strength of the hopping between sites k and I. 
The creation a k and annihilation hk operators follow the standard commuting rules for 
a discrete model given by equation (jl.30|) . 

The Hamiltonian above is quadratic and can be decoupled in order to get the ground 
state. Additionally, the Hciscnberg equations of motion for the creation operators pro- 
duce a complete set of differential equations that can be written as, 



N 



N 



N 




(3.1) 
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dt 
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or cquivalcntly, 





-Rl,2 


-Rl,3 


-K 1,2 


R2.2 


R2.3 


-K 1,3 


D* 

-n- 2,3 


&3,3 



R* 



1,N 



R' 



2,N R*3,N 



R* 



N-1,N 



Rl,N 




( 


01 


R2,N 






&2 


R3,N 






&3 


Rn-i,n 








Rn,n 


) 


V 





, (3.2) 



d a 

dt 



iRa, 



(3.3) 



where oik describes the corresponding creation operator in the Heisenberg picture, 



a k = e- itH ale itH . 



These equations are complemented by the initial conditions, 



(3.4) 



a k {t = Q) = a\. 

The state, on the other hand, is given in terms of the Heisenberg operators by, 



1 N 



1 1 rijl k 



(3.5) 



(3-6) 



The way in which particles are distributed across the chain is determined by the con- 
stants rij. 

From equation (|3.3p it is in fact possible to obtain the ground state of the system through 
imaginary time analysis (equation (|4.7|l ). In doing so, the ground state is found to be, 



A' 



M 



\g)= E c <4 i°>< 



(3-7) 



\k=0 



where coefficients Ck are the components (possibly complex) of the ground eigenvector 
of matrix R, so that, 
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N 

£M a = 1 - (3-8) 

fc=l 

From equation (|3.7[) we can extract some information. Nevertheless, because the de- 
scription of the state demands exponentially growing resources which scale with both TV 
and M, calculations involving non-local correlations can be fairly challenging. In order 
to set state (|3.7[) in a way that can be easily handled, we will apply unitary operations 
to the state so as to simplify it as much as possible. The idea is to reduce \G) to an 
expression that can be written using MPS. Then, if the unitary operations only involve 
transformation to first neighbours, the updating method presented in previous sections 
can be applied. This would provide us with a complete description of the state that we 
can efficiently use. 

We first operate locally on individual sites using the transformation, 

In this equation 6k is the phase of the complex number c^.. As a result of this transfor- 
mation we get, 

4 H- e-**a£. (3.10) 

In such a way that the complex phase is cancelled out. Additionally, any operator 
different from &k remains unaffected by the transformation. Once this has been done 
for every operator &k, the new coefficients Cfc in equation f|3.2[) are real. 
Subsequently, we apply unitary operations involving first neighbours in the following 
fashion, 

(j( 2 ) = e ~ i <t>k(ii{a.[ + l a k -ala k+1 )) _ e ~i<t>kQk ^ (3-11) 

In physical terms, operator Qk is known as the current, and its eigenvalues indicate how 
many bosons circulate in between sites k and k + 1. As a result of this unitary operation, 
the pair {aj^etl+i} transforms as, 
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c k+ ia\ +l + c k a\ -> (3-12) 

h 

2 



c k+1 cos I J - c fe sin ( ) ) oL ! + I c fc+ i sin I I + c fc cos 



Therefore, we can cancel operator a\ +1 by just choosing an angle 4> k satisfying, 



(3.13) 



Consequently, in order to reduce the state to one single mode operating on the vacuum, 
we apply the operations (|3.11|> . starting from k = N — 1, to every couple of consecutive 
operators. Note that every time a transformation acts on any pair of modes, the coef- 
ficient that accompanies the creation operator that is not taken out is affected, but it 
always remains real. 

When the reduction is completed, we are left with a state as (up to an overall constant), 



Iff) 



1 



M 



|o>, 



(3.14) 



which can be easily written using MPS. The next step consists in applying the inverse 
operations in reverse order to the state using the method presented in section ll~3l From 
now on, we will refer to this reduction operation as state folding. 

There is also the case when we must deal with numerous summations acting on the 
vacuum. For instance, in chains with an initial state given by bosons arranged on 
different positions. In this kind of situation state folding can be utilized too. Let us 
consider the case when just two summations get involved so that the state reads, 





M 2 




Mi 


N 




N 










|0>, 


*;=i 




1=1 




V J 




{ ) 





(3.15) 



where M\ is the number of bosons originally allocated on one site of the chain and Mi 
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has an analogous meaning. We can use the standard technique to fold S±, but it is worth 
saying that in the process the coefficients of S2 are also affected. Next, we can fold the 
new 5*2 from ajy just until a 2 , since folding a\ in S2 would unfold Si. As a result, the 
folded state must be written making explicit reference to such last folding operation in 
the form, 

\g) = (a\) e-^(a\) |0>. (3.16) 

Consequently, \g) can be written in MPS by first translating \Mi, 0, 0) into MPS and 
then applying e -1 ^ 1 ^ 1 . In this way, we are left with a canonical decomposition that 
can be used to apply the last operation in the row. With this in mind we set down the 
state with explicit reference to the local coordinates of the first position and then apply 
[a[j in a very straightforward way. In doing so we get, 

l.9> = E r i?$ ] V(3 + M 2 )l\j + M 2 ) | 7 ) ^ . (3.17) 

7=1 

In latter equation we have made use of the fact that the number of bosons in the 
local basis of the first site \j + M 2 ) is determined by the number of bosons in the 
complementary Schmidt vector I7) and therefore there is only one relevant coordinate 
that describes the local basis. This explains the lack of a summation symbol for the 
label j. From the expression above the canonical coefficients of state \g) can be directly 
obtained, namely, 



A 7 = Vti + M 2 )l\ i y + normalization, (3.18) 
r ,[i]j+M 2 = r Mi_ (31Q) 

However, because this is not a unitary transformation, it is important to show that the 
canonical decomposition obtained in this way is consistent. It is not difficult to see that 
the canonical tensors attached from the third site onwards do not suffer any modification 
whatsoever since this section of the chain is made of complementary partitions which 
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Expectation Value of the Number of Particles 
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Figure 3.1: Expectation value of the number of particles. Top plots are grey maps of 
the lower plots. Dynamics is generated by turning a potential barrier in the middle on 
(left) and off (right). 



contain only one Schmidt vector. This is because there is no boson at all between the 
third and the last position. Wc can see that such is indeed the case by noticing that the 
operations performed on the vacuum involved only the first and second places, thereby 
only these two positions can hold any non-vanishing population. On the other hand, we 
know that the tensor elements in positions one and two depend directly on the Schmidt 
vectors {|7i)^}- These vectors all have well defined quantum numbers on account of 
particle conservation. Therefore, when we apply I a[ ) we are actually lifting the boson 
occupation which means that the resulting vectors are valid Schmidt vectors which are 
orthogonal to each other. 
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3.1 Quench in trapped systems 



We now show how the method presented in the previous section is capable of delivering 
results in challenging scenarios. In what follows, we assume that energy is given in terms 
of the recoil energy, 

hk 2 

Er = ^' (3 - 20) 

where k is the wavelength of the confining laser and m is the atomic mass. Let us 
consider a chain of 100 sites and 100 bosons which has been cooled down to the ground 
state in the presence of a trapping potential given by, 



R iA = Q(i - 50) 2 , (3.21) 

with ft = 0.00046i?R, a reasonable experimental factor according to reference [5^. Ad- 
ditionally, we assume that particle exchange can take place among any two places in 
the chain, not only next neighbours. The intensity of the hopping is proportional to the 
off-diagonal matrix elements, 

Rl ' J = \i-jV < ' 3 ' 22 ' ) 

for every i ^ j and with 3 = 0.3-Er. This choice of hopping can be justified in chains 
with long-range exchange effects, such as expected in situations where the Coulomb 
potential plays an important role. In a first set of simulations, we find the ground vector 
of matrix R and insert the coefficients in equation p.7[) . Dynamics is generated by 
instantaneously turning on (off) a very high potential barrier in the middle of the chain. 
This barrier is written as, 



Ri,i = 1000£^ for 45 < i < 55. (3.23) 

Next, we find the ground state of the system including the potential barrier so as to 
generate a distribution with two condensates at the sides of the wall and then turn 
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the potential barrier off. Figure 13.11 shows the behaviour of the expectation value of the 
number of bosons for the above mentioned configurations. As a result of the quench very 
characteristic wave patterns are generated depending on how the quench takes place. In 
both cases, however, fleeing waves are generated after a particle bulk assembles in the 
chain centre. In the first case, when the barrier is turned on before the condensate is 
released, we can see that a large part of the condensate remains in the middle of the 
chain, just where the barrier stands. The barrier domain seems to determine a zone 
from which bosons can hardly escape. It is reasonable as the high potential difference 
between sites inside and outside the wall zone energetically hinders any quantum jump. 
In contrast, the pattern generated in the case when the dynamics is generated by switch- 
ing the potential barrier off is more compatible with uniform expansion. Nevertheless, 
in both cases we can appreciate the effect of the inclusion of long range hopping terms 
in the Hamiltonian. When the bosonic waves travel outwards they gradually lose par- 
ticles. However, this loss does not translate in dissipation, instead, a new boson bulk 
reassembles around the centre of the chain and a new expansion starts again. This effect 
is clearer in the case of expansion without intermediate potential, but it also takes place 
in the other case analysed. Such profile is in contrast with the behaviour of travelling 
packets of chains with next-neighbour hopping only. In the latter case the changes in 
the form of the wave pattern develop locally around the boson bulk, and the effect of 
particle exchange between distant places of the chain is much less pronounced. 



3.2 Entanglement as a result of collision 

In this section we want to show how the two-sum folding method can be used in time 
dependent problems. Let us suppose that in a boson chain we have two very well 
localized boson clouds at the ends. We then use a physical mechanism to accelerate the 
clouds and make them interact with each other. After this collision process, the scattered 
bosons are collected and taken back to the chain terminals. Because of the interaction 
between the boson packets, the collected particles at the ends are strongly entangled, as 
discussed in [SHI HOI ED E2 ■ Certainty, this entanglement would be potentially useful for 
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Figure 3.2: Sketch of a bosonic collision 

multiple quantum information procedures. Before focusing on entanglement, however, 
we will look at how the bosons get transmitted from one end of the chain to the other 
in the presence of a perturbation in the central part. In such a system the Hamiltonian 
would be ideally given by equation (|3.1[) with a matrix R that can be written as, 

R= j x + ee-P J >, (3.24) 

where J x ,y,z are the standard angular momentum operators, e, on the other hand, rep- 
resents the intensity of a symmetrically localized perturbation in the middle of the chain 
and /3 measures the spread of such perturbation. One reason to put matrix R in terms 
of the angular momentum operators is that in this way we can highlight the analogy 
between perfect transmission in a boson chain and the physics of angular momentum. 
In fact, in the absence of a term proportional to e, this would be basically a particle 
spinning around the x-axis. The perturbation has been deliberately chosen in the form 
of an exponential in order to generate a decaying profile as going from the middle of 
the chain towards the terminals. This form is also very convenient to perform analytical 
calculations as we will show. In addition to this formulation of matrix i?, we must estab- 
lish a relation among the state kets in the Hilbert space and the actual operators. This 
can be worked out by just comparing the effect of the angular momentum on the corre- 
sponding kets with the predicted behaviour of the operators in the perfect transmission 
scenario. From this we can infer, 
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\j,m) ^ a]_ m+v (3.25) 

In this expression we used the standard notation for the angular momentum kets, that 
is, \j, to) are the eigenstates of the z-componcnt operator J z . The quantum number j is 
related to the total number of sites in the chain by the identity, 



i = ~TT~- (3-26) 



In this way, the evolution of the Heisenberg operators can be studied by following the 
dynamics of the associated kets. In a problem where the physics of the system is dictated 
by equation p.24[) the wave function acquires the following form after half a period 
evolution, 

m = n)) =e-™ (J * +ee ~''%,j). (3.27) 

Here we have implicitly assumed that the length of the angular momentum j is an 
arbitrary integer which accounts for N being odd. We have focused on the evolution of 
\j,j) because that is the ket relevant to operator a\ and therefore the one which contains 
information about the transmission of bosons initially prepared on the first site of the 
chain. Making use of time dependent perturbation theory we can expand the latter 
expression in the following way, 



= 7r)«e-* r %j>-ie|e> ) 
\e) = ft dte-^-QS'e-rte-* 3 '^, j), (3.28) 

of course, this approximation holds as long as the amplitude of the dynamics generated 
by the integral in the second term remains small compared to the unperturbed evolution. 
For this to happen it is important not only that e is a little fraction but also that the 
dynamics develops out of resonance. Nevertheless, in this particular case where we 
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restrict ourselves to well defined time intervals, the emergence of cooperative resonances 
is quite improbable. So, the integral above can be manipulated as, 



|e)= / dte-^-^e-^e^-^e-^UJ). (3.29) 
Jo 

From the properties of the angular momentum we know 



-17T J, 



\j,j) = \j,-j)- (3.30) 



Further simplifications apply by noticing that the three first exponentials from left to 
right underline a unitary transformation. As a result we can write, 

|e)= / dte-^^^-^-iy^lj,-:)). (3.31) 



JO 

Because the operator on the integral is squared, we do not have a way to translate the 
expression into a sum of kets. Therefore, we rather use the well known result from 
integral calculus, 



oo 

7T b 2 



dxe _ {ax * +bx) = l e is (3.32) 
V a 

to put |e) in a more operational fashion, 

dxe~i ( d te «(^(t)4-«n(4)i„)|^_^_ (333) 



o 



In this way, we can now work out how the exponential operator acts on the state by 
making use of Schwinger's alternative representation of the angular momentum |73j . 
For this we define a couple of independent sets of bosonic operators (contrary to the 
operators that describe the bosons on the chain, these do not have a direct physical 
correspondence) , 
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= 1, 



a_, 



= 1, 



a + ,al =0, a_,at_ =0 



(3.34) 



It can be shown that these operators form a structure that correctly reproduces the 
angular momentum through the following equivalence transformations, 



along with, 



a\a,-. 



J- = a_a+, 



J z = ^(a + a+ — al<z_), (3.35) 



Cot \j+m(fA )j-m 

\j,m) = ^= == |0+,0->, (3.36) 
V(j + m)!(j _ m ) ! 

for the state kets. In terms of the bosonic representation the perturbed state reads, 



OO 2 



^47r/3(27r)! 

dfe «(£2iW( a |a + -ala_)-i^(a|a_-aL5 + )) _^ 2j | 0+j 0_). (3.37) 



We now group terms in the second integral as indicated below, 



(fit) ' |0+,0_) = (e^ale-^) ' |0 + ,0_), (3.38) 
where, 



A = 
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and 4> is just an auxiliary parameter which can be set to (j> = 1 whenever is conve- 
nient. This identity could be verified, for instance, just expanding the right hand side 
of equation Q3.38P and noticing that, 



- a *|0+,0_) = |0 +) 0_>, 



(3.40) 



which can be seen if the exponential on the left is expanded in power series. It can 
be seen, therefore, that finding an analytical expression for the evolved operator is 
structurally analogous to the formulation of the free bosonic model of section[3J In fact, 
here we also get a system of equations that looks like, 



a+ 



a- 



( 



— f cos(i) ifsin(t) 



-i% sin(i) 



cos(i) 



a,- 



(3.41) 



along with the initial condition, 



"+(0 = 0) = a+, 

&_(<£ = 0) = a_, (3.42) 

and the usual definition Therefore the complete expression for the 

evolved operator is, 



d_ = e'i ™M fsin ( & l + cos f^R) a t ) 



Inserting this result in equation p.37[) we obtain, 



/UU .22 r " 

dye' 1 ^- / dte ij2xcoB ^ (3.44) 
-oo Jo 

»t . /0'2/sin(i)\ , f / jy sin(t) \ \ 2j 

<sm +aLcos |0+,0_), 



where we have performed the variable change x = jy. Expanding the binomial and 
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applying the boson operators on their respective vacuum states we get, 



^ fc=0 v v 



' 72/ sin(i) \ 2j fc / jy sin(i) \ fc . 
sin cos |2j-fc + ,fc_>. (3.45) 



Next, we approximate the trigonometric functions by the first term in their series ex- 
pansion, 



sm 



jysin(t) \ jysin(t) 
2 2 ' 



^ jy sin(t) ^ 



(3.46) 



This approximation is valid as long as the negative quadratic exponential in the first 
integral of equation (|3.45[) suppresses the contribution of the trigonometric functions for 
large arguments. Such is indeed the case for a wide range of values of /3 and j, since 
the argument of the first exponential decreases quadratically while the argument of the 
second exponential scales linearly. After some algebraic manipulations we arrive to, 



j dte ij2ycos{t) sin(t) 2 ^ fc J |2j - k + ,k-). (3.47) 

The integral in curly brackets can be worked out in terms of Gamma and Bessel func- 
tions. So we are left with, 



W = ffi ( V + 5) £ (I) ^ - ^ 

(3.48) 
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We now use the first term in the series expansion of the Bessel functions around x = 0, 



which provides a good estimation in the range of small arguments. We expect this to be 
a good approximation since the negative exponential in (|3.48p attenuates contributions 
from large values of the argument. The integral is therefore reduced to, 

(3.50) 

From the integral above wc can infer that to first order the contributions corresponding to 
odd positions all vanish. Consequently, through direct integration and some re-indexing 
we obtain, 

w = /§) (If 121 - 2 ' - im - 2,+> ^ (3 ' 51) 

Additionally, we utilize the following identity, 

r(n+0 =^(2n-l)H, (3.52) 
which helps simplify the whole expression to, 

(3 - 53) 

Now we can use equations (|3.36[) and (|3.25[) along with the original perturbativc expan- 
sion in equation (|3.28D to establish a relation for the bosonic operators, 
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Figure 3.3: Behaviour of function F(q) for j ' = 5 and j3 = 5.5(arbitrary units). 



In order to illustrate our results, we have depicted in figure [3731 the continuous function 
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r(2j + 1) 



r (;-* + !' 



r(2j-2 g +l)r(2g+l) T(j-q + l) 



(3.55) 



which underline the discrete function in the sum of equation (|3.54l) . Equation (|3.54l) 
shows that when j3 = the evolution is determined by the integrable part of (|3.27[) 
and therefore only the creation operator a N survives. This means that particles are 
efficiently transferred across the chain terminals. For intermediate values of (3, on the 
other hand, the distribution function F(q) gets delocalizcd and particles spread across 
the chain. The most interesting case, however, occurs for large values of (3. Such as it 
is shown in figure 13.31 the function is highly localized around values of q corresponding 
to a\. As a result, bosons will not be spread all over the chain any more, instead, only 
the terminals will be macroscopically occupied. As this happens only for large values of 
/?, we conclude that this is characteristic of chains with tightly localized perturbations 
around the centre. One can compare this situation with the dynamics of an individual 
particle in one dimension in the presence of a potential barrier. From this parallel we 
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Figure 3.4: Terminal occupation after half a period time. On the left we fixed the number 
of bosons and vary the length of the chain. On the right we maintain the number of sites 
fixed. In this situation the curves corresponding to different number of bosons collapse 
into one single characteristic curve. 



can think that when the particles get over the middle of the chain the perturbation 
simply acts as a thin wall that causes reflection without altering the wave packet shape. 
As a result, reflected particles preserve the coherence necessary to drift back into their 
original chain terminal while transmitted particles follow their predetermined path. This 
special feature is of potential usefulness in scenarios in which particle localization plays 
a crucial role, such as the one we now focus on. 

Let us now consider a situation in which particles are originally prepared in a separable 
state with all the bosons distributed evenly between both chain terminals. We use the 
PTH profile to induce coherent particle transmission but we also take into account the 
interaction of the boson clouds around the centre of the chain (figure 13. 2|) . Formally, 
the evolution is dictated by equation (|3.1[) with the following coefficients, 



1 



Ri ji+l = -y/i(N-i) 



(3.56) 



and, 



#f,f = = M- (3.57) 

Here we assume that the number of sites across the chain is even, since our algorithms are 
designed for such particular configurations. Hence, we assume that the interaction takes 
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place in two central sites, but in a chain with N odd, a highly localized perturbation 
can be modelled in one single central place. Recall also that the previous analysis 
corresponds to chains of odd size, but we expect that our results are robust against 
this small discrepancy since the process we study is quite general and straightforward. 
Evolution is simulated applying the two-sum folding method to the state, 



\ip(t = tt)) = &i(t = ir)-& N (t = tt)-|0), (3.58) 

where operators oti and come from solving equation (|3.ip for the coefficients in (|3.56l) 
and (|3.57[) . After the state is written in MPS, we find the reduced density matrix of the 
ends of the chain. This matrix must be written in the basis of the number of particles 
in order to get the logarithmic negativity associated with the state (sections [T] and [5]). 

We first observe how efficient the collection process is. Figure [33] depicts the proportion 
of particles on the ends once the process has taken place, namely, 

(a[ai) + (ajy-ajv) 

(3.59) 

From both plots we can see that the particle proportion is very close to 1 which means 
that almost all the bosons are collected on the chain ends along the domain of [i. This 
corroborates the analytical results previously found. In figure 13.41 we use the derived 
variable ^ because it appropriately rescales the plots in suitable proportions and puts 
all the results in a similar scale. This allows us to see that the collection efficiency 
notably improves when the size of the chain goes up, which is reasonable as in this way 
the perturbation gets more localized in relation to the rest of system. For a chain of 100 
places, collection efficiency remains above 99% independently on the amount of bosons 
involved. This would be a very positive aspect in a scattering experiment, where it is 
important to have well localized boson clouds in order to measure local observables. 

Once the bosons settled down in the chain terminals, the system is very entangled as 
a result of the collision. We can understand this entanglement as coming mainly from 
the quantum correlations among the dynamical degrees of freedom at the moment of 
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Figure 3.5: End-to-end entanglement (EEE) of a chain of 100 sites for t = it. 

the interaction, as it is in the centre of the chain that the bosons become more speedy. 
However, such entanglement is difficult to use straight away, not only because the particle 
distribution is very diffuse, but also because it is difficult to distinguish between highly 
overlapping boson clouds. Therefore, it is better to wait until the system comes back to a 
zero momentum state, where two static boson bulks lie well separated from one another 
and entanglement becomes accessible. From figure [3~5l we can see that the entanglement 
generated by this mean is quite substantial for optimal parameters. Additionally, it 
steadily grows according to the number of particles involved in the collision. When 
the interaction parameter is fi = 0, the boson packets just get across each other and 
no quantum exchange takes place. Under this circumstance the possibilities of any 
entanglement emergence are null. Nevertheless, as fi is slowly turned up, entanglement 
grows quickly. In this case the boson clouds interact gently, but this small perturbation 
over the integrable path effectively scrambles the degrees of freedom and entangle the 
boson waves. Notice that the entangling process takes place mainly in the middle of 
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Figure 3.6: Logarithmic dependence of the maxima of figure [ 



the chain, when the packets are half the way between the chain terminals. Once the 
interaction is over, however, the particles will travel towards one of the chain corners, 
depending on whether they were transmitted or reflected in the collision. Interestingly, 
due to the absence of interaction the entanglement between the boson bulks does not 
build up during the travelling, instead, it is kept constant. 

On the right side of the graph in figure 13.51 we can observe how entanglement decays 
for large values of the interaction constant. Logically, when the interaction is infinity, 
bosons just bounce back in the middle of the chain and return to the original terminal. 
As a result, there is no way for places in opposite sides of the chain to "communicate" 
and neither correlations nor entanglement do arise. However, the entanglement decay 
displayed by the system for large [i is much slower than the pace at which entanglement 
soars close to \i = 0. 

Figure |3~51 shows that the best entanglement production takes place in the range 0.25 < 
j7 < 0.75. Additionally, the positions of the maxima are confined to a small interval, 
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making clear that the process is to some extent generic. Likewise, the values of the 
maxima keep a logarithmic relation to the number of bosons, such as can be appreciated 
in figure 151)1 Numerical fitting yields, 



Emax = 0.65 + 0.79log 2 (M). (3.60) 

This basic behaviour can be explained by formulating the problem in a slightly simplified 
way. As the entangling mechanisms are highly localized in the middle of the chain, we 
can assume that the amount of entanglement generated at the moment of the collision is 
equivalent to the entanglement between the ends when the process finishes. In a chain 
without interaction, the state would be given by, 



|V> (<=£)> = *!(*=!) 2 2 |0)=a A »= £<*4 |0>- (3.61) 

\fe=o / 

This result comes from the fact that the particle distribution is completely symmetric 
for a time equal to a quarter period. Consequently, in a chain with no static interaction 
the Hcisenbcrg operators became identical and the state can be written as a single sum 
acting on the vacuum state. As the focus of the analysis is entanglement, we can operate 
on the state using unitary transformations that preserve the amount of entanglement 
contained in the system. Additionally, because most of the particles end up on the ends 
of the chain and also because entanglement arises mainly on the interacting zone, we 
can assume that the amount of entanglement between the halves of the chain represents 
a good estimation of the entanglement between the terminals at the end of the process. 
Therefore, we apply the one-sum folding technique introduced in section but this time 
we do not reduce the expression to one single creation mode. Instead, we operate on 
pairs of operators belonging to different chain halves, starting with the pairs of operators 
on the ends and going towards the operators in the centre of the system. The result of 
this reduction can be written simply as, 
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Norm x (a f „ + a f „ , . ) M |0) . (3.62) 

2 2 + 1 

Entanglcment between the chain halves results from expanding the binomial, 



Z\lww^W» lM - km ' (3 - 63) 



k=0 

which provides us with the Schmidt coefficients of a symmetric partition, 



/ 1 M! / 2 \ 4 i /, «v , 

so that the coefficients can be approximated by a normal Gaussian. In this auxiliary 
model Log- negativity is given by |62j , 



En = 21og 2 (j> fe ) « 21 °S2 ((^) 1 / da;e ~^ ■ ( 3 - 65 ) 
Finally, after some direct calculations we obtain, 



E N w log 2 (27r) + log 2 (M) = 2.65 + log 2 (M). (3.66) 

Which positively verify the logarithmic behaviour of entanglement previously underlined 
by the numerical analysis. The discrepancy in the numerical factors comes as a result 
of the approximations involved in the derivation. 

In summary, we have shown how our folding technique can be used to simulate the 
ground state as well as the dynamics of highly entangled systems in an efficient way. 
Here we chose to perform computations in problems of physical relevance such as quench 
and scattering but the method can be applied to a wide range of problems. In the case of 
quench, for example, we found that the wave patterns predicted by the simulations are 
characterized by the propapagation of bosonic waves all over the chain. For the case of 
scattering, on the other hand, we showed how boson packets get strongly entangled after 
they collide and are then collected back on the chain ends. As the maximum amount 
of entanglement grows logarithmically with the number of particles in the system, the 
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configuration proposed has potential application in experiments in optical lattices with 
direct implications in quantum information processing. 
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Chapter 4 

The Bose-Hubbard model and 
time evolving block decimation 



The BH model describes a collection of interacting bosons that can hop among neigh- 
bouring sites and undergo on-site repulsion when more than one boson occupies the 
same site simultaneously. The BH model in one dimension is given by the Hamiltonian, 

N jj N-l N 

H = — !) — Jfc ( a fc+i afc + ata k+1 ) + ^2/i k ala k . (4.1) 

fe=i fc=i fc=i 

As usual, the bosonic operators obey the commuting rules introduced before in equation 
(|1.30[) . while N is number of sites in the chain and the implicit integer M is the number 
of bosons. The set of Uk determine the intensity of the repulsion and the set of Jf. 
determine the intensity of the hopping across the chain. Finally, the set of /ifc account for 
the chemical potential. The BH model has been extensively studied since the discovery 
of superfluity in 4 He |J. Although exact analytical solutions are not available, a lot 
of insight can be obtained out of careful inspection and clever mean field approaches. 
One important aspect of the study of the BH model is that such study is frequently 
carried out following a grand-canonical ensemble approach. This means it is assumed 
that the boson chain is being constantly fed by a powerful supplier that continuously 
provides bosons as well as energy in order to maintain the thermodynamic equilibrium. 
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Figure 4.1: Phase diagram of the BH model. Inside the Mott insulator lobes the number 
of particles per site remains fixed. 

Nevertheless, it is important to mention that here we do not follow such description. 
Instead, our chain is fully isolated from any outside disturbance and the system state is 
given by the ground state of the Hamiltonian or by real-time unitary evolution. In figure 
14.11 the phase diagram of the BH model interacting with a supplier bath is presented. 
As can be seen, the model is characterized by a supcrfluid phase and a Mott insulator 
phase. The difference between the phases is quite clear. In the superfluid state, bosons 
display non-vanishing hopping-scopc-length. We can say that this phase is driven by 
hopping, since J » U . On the other hand, in the Mott insulator state bosons remain 
tightly fixed on their hosting positions. In this phase hopping is highly hindered due to 
the loss of energy balance when the Mott insulator configuration is slightly disturbed, 
which results from U >> J. Formally, the Mott insulator phase is characterized by the 
following properties, 

• Integer or commensurate density. 

• Existence of an energy gap. 
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• Zero compressibility. 

• Exponentially decaying correlations (ajcij). 

The first property is clear from the phase diagram of figure 14. f I The number of bosons 
per site is fixed inside the lobes that enclose the Mott phase. Commensurate densities 
correspond to configurations where the number of bosons is not an integer multiple 
of the number of sites, but a rational multiple. Over the y-axis and in between the 
lobes, the Hubbard model describes a superfluid, since the total number of particles 
fluctuates between two fixed values. The second property is of outstanding experimental 
importance, since the energy gap allows us to identify the Mott phase in the laboratory 
by spectroscopically probing the sample [5]- Additionally, the gap can be used in EIT 
experiments, where very long storage times have been observed [SJ. The gap refers to 
the energy cost associated with making the bosons in the Mott state to jump out of 
their original positions. The third property means that the density of the system does 
not change as the chemical potential is varied, or equivalently, = 0. This also can 
be inferred from figure 14.11 by noticing that if we move from inside one of the Mott 
lobes in the direction of the vertical axis the number of particles would remain fixed as 
long as we stay in the Mott insulator state. We can think of this as a system in which 
the number of particles remains fixed despite the channels that define particle exchange 
between the system and its environment remain open. 

Conversely, the superfluid is characterized by particle tunnelling all across the chain. 
The main properties of the superfluid are, 

• Can exist at any filling. 

• Eminently gaplcss. 

• Finite compressibility. 

• In one dimension, power law decay of correlations (ajaj). 

The first three properties stand opposite to the properties of the Mott insulator. The 
fourth property establishes that the superfluid is liable to develop correlations at long 
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scale. This can be understood physically by noticing that as bosons jump from one place 
to the other, they distribute information all over the system. This property constitucs 
the basis of the schemes that will be proposed in subsequent sections in this work. Both 
the Mott insulator state as well as the transition to a supcrfluid have been verified using 
cold atoms in optical lattices [51 1141 17]. In a typical experiment, atoms are taken into 
a magnetic trap and then cooled down using an optical lattice of rctroreflcctcd diode 
lasers. This creates an arrangement of atoms where the resulting optical potential depths 
Vx,y,z are proportional to the laser intensities and can be expressed in terms of the recoil 
energy, 

hk 2 

E B = ^, (4.2) 

with m the atomic mass and k the wavelength number. To prepare ID arrays, two lattice 
lasers are given high intensities in such a way that hopping can only efficiently take 
place across one axis [7J. In terms of experimental parameters, hopping and repulsion 
coefficients arc given by, 



(4.3) 



and, 



d V En \E R / 

where Vq is the axial lattice depth, V± the depth of the lattice in the transverse direc- 
tions, a s the s-wave scattering length, d the lattice spacing. Capital letters represent 
experimental constants established by the geometry of the system. For example, for the 
experiment reported in reference [M] such constants are [TH [39] , 



A = 1.397, B = 1.051, C= 2.121. (4.5) 

Spatial variations in U and J can also be implemented using detuned lasers sent through 
specific sections of the lattice as shown in reference [TT] . 
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Due to the vast number of quantum states necessary to span the Hilbert space, we must 
utilize MPS in order to be able to handle chain sizes of physical relevance. In an isolated 
boson chain, the conservation of the total number of bosons along with bosonic statistics 
determine the size of the basis needed to write an arbitrary state, namely, 

basis size = ^^ r ,"!" ^ — rr-. (4.6) 
M\(N-l)\ y ' 

In practical terms, the amount of resources needed to simulate the system grows expo- 
nentially with the number of sites and bosons. Studies in boson chains using conventional 
bases have been done in chains of up to 10 sites and 10 bosons [75] . 
In order to find the ground state of the chain, we apply imaginary time evolution to a 
given initial state written in MPS. Such an initial state must have some overlap with 
the system ground state. The imaginary evolution is explicitly given by the formula, 

\^ G )=hm T ^ 00 „ e ~^°|„ . (4.7) 

Additionally, we use second order Trotter expansion, which consists in splitting the 
evolution operator in a product of non-commuting operators as, 



Each operator in the above product can be split in a product of several commuting 
operators, each corresponding to a specific pair of neighbour sites as shown in figure 
14.21 After every updating step, in which the state is evolved using the Trotter formula, 
we retain all the Schmidt coefficients greater than A^ x 10 -14 , where A J is the greatest 
coefficient at site k. Then we use the Schmidt coefficients to find out the canonical 
representation. The maximum number of Schmidt coefficients in the chain is denoted 
by the letter x- The larger \ the more expensive it is for us to simulate the chain. In 
our program we use time slices in the range 10 -5 < St < 10~ 3 depending on the chain 
size. Usually, simulations in long chains as well as very entangled chains require the 
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Figure 4.2: Sketch of operator splitting. 



smallest St. Ground state convergence is evaluated through the criterion, 

|1 - (V>(t)|V>(t + 5t)))\ < 1CT 14 . (4.9) 

The accuracy of the final state depends on factors such as the length of the time slice, the 
original state, and the parameters in the Hamiltonian. In figure l4~4l we show simulations 
for trapped bosons in chains with different repulsion. The chosen regimes are compatible 
with experimental observations of Mott insulator (driven by repulsion) and superfluid 
(driven by hopping). The original state is a separable state for which the canonical 
decomposition can be easily written, namely, a state with one boson in each site. In this 
case the canonical decomposition can be written as, 



A'™ 1 = l,Vn 

if;; 1 = #,vn. (4.10) 

For this state x = 1- In a standard simulation, the maximum number of Schmidt 
coefficients \ gradually steps up as in figure [CU In most cases this variable saturates, 
but sometimes it instead peaks and dips, then reaches an equilibrium value. This is 
particularly the case for PTH chains with no repulsion. Figure 14.31 shows the behaviour 
of x m imaginary time evolutions corresponding to different time slices. From such graph 
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we can see that simulations with large St require large X- However, the number of steps 
necessary to reach convergence is fairly low. On the other hand, when the size of St is 
reduced, the maximum \ goes down while the number of steps necessary to converge 
the state goes up. We can actually establish the appropriate value of x by comparing 
the energy of the obtained ground state with the exact analytical results. We know that 
the energy of a PTH chain with M bosons is given by, 

Eg = -A (^Y 1 ) M. (4.11) 

Similarly, we can estimate the energy Eq of the state written in tensor notation from 
the change in the state norm assuming that the state itself is not hugely different from 
the ground state, 



e-^c) =e- 5TE °\4, G ). (4.12) 

Table |4~T1 shows the energy values obtained in this way. We can see that states obtained 
using relatively big St slightly underestimate the ground energy. In this case we can say 
that the theoretical error involved in decomposing the evolution operator in a product 
of small unitary transformations is more influential than roundoff errors. Additionally, 
we can see there is an optimal value of St for which the ground state energy contains the 
maximum number of significant correct figures without it going below the actual ground 
energy. Similar analyses in chains of different sizes show that in PTH chains optimal 
convergence is achieved using St = 10 N . In the general case, the optimal size of St is 
inversely proportional to the scope of the boson tunnelling in the ground state, which 
depends on the intensity of the hopping. Hence, for arbitrary hopping and repulsion 
profiles we can write, 

10~ 3 

St = ^-, (4.13) 

where N e is the effective tunnelling length of bosons in the state, which must be esti- 
mated out of physical insight. 
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Figure 4.3: Ground state convergence in PTH-repulsionless chains for 10 sites and 10 
bosons. These simulations correspond to the energy values in table |4"~1"1 



St 


E G 


0.1 


-.902854035422649E+02 


0.05 


-.900740540911757E+02 


0.01 


-.900029984667920E+02 


0.005 


-.900007499039339E+02 


0.001 


-.900000299962281E+02 


0.0005 


-.900000074855496E+02 


0.0001 


-.899999999159858E+02 


0.00005 


-.899999988644993E+02 



Table 4.1: Ground state energy in PTH-rcpulsionless chains with no repulsion for 10 
sites and 10 bosons. These values correspond to the energy graphs in figure [431 The 
value closest to 90, the actual energy, is written in boldface. 
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Figure 4.4: Simulations in boson chains of 100 sites and 100 bosons using TEBD. The 
error bars in the lower graphs represent fluctuations. Left: Uk = 2, ~ 0.14, fj, k = 
0.0046(fc - 50) 2 and St = 0.00001. Right: U k = 0.42, J k = 0.14, /x fc = 0.0046(fc - 50) 2 
and St = 0.00005. All these constants are given in arbitrary units. 



79 



As it has been pointed out in several studies [751 E3] > it is the variable % that ultimately 
determines the simulation efficiency. This fact is quite reasonable, on the one hand \ 
defines how much memory is necessary to store the canonical decomposition, and on the 
other hand it determines the number of numerical operations that must be performed 
in order to upgrade the state [221 163] . Similarly, this number is directly linked to the 
entanglement profile displayed by the state. As can be seen in figure POI states showing 
more entanglement also develop higher \. Therefore, practical simulations can be done 
in systems where analytical correlations are not too strong. In the BH model particu- 
larly, we know that close to the Mott insulator regime quantum correlations are quite 
moderate. Additionally, the superfluid state is characterized by certain regularity which 
somehow keeps x from growing indefinitely. However, it is the state in between the 
phases what brings more challenge. In this case, a kind of strongly enhanced entangle- 
ment embraces the state and x grows well above manageable limits. According to mean 
field theory, the transition from Mott insulator to superfluid in the thermodynamic limit 
takes place at y = 11.6 [H[S] for a chain with unit filling. Figure POl shows simulations 
on each side of such transition point. These graphs show the characteristics of the Mott 
insulator and the superfluid. As can be seen, the expectation value of the number of 
particles in the Mott state shows integer density over a wide range of positions as well 
as small fluctuations. The superfluid features large fluctuations, which anticipates the 
existence of long range correlations. In general, the simulations of the superfluid require 
more computational resources than the Mott state. As discussed before, the Mott state 
is characterized by the existence of a gap as well as integer density. The superfluid on 
the other hand is associated with intense boson hopping and long range correlations. 

Similarly to ground state simulations, the performance of TEBD in the real time con- 
text is highly dependent on the variable x- However, the development of this variable 
in dynamical problems is different from the one observed in imaginary-time evolution. 
While for the latter case x saturates to a reasonable value, for the real time case it grows 
steady, the same that the number of resources necessary to describe the state [77]. As 
a consequence, numerical dynamics using TEBD is usually restricted to short times 
regimes, where the size of the canonical decomposition can still be handled. Alterna- 
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Figure 4.5: Dynamics of a boson chain using TEBD. The system is prepared in a Mott 
insulator state with one boson per site. The chain parameters are characteristic of the 
supcrfluid phase, namely, Uk = OAEr, Jk = Er and = 0. In order to get reliable 
results, we use a very thin time slice St = 0.00005Eji/h. As a way of keeping the size of 
the canonical decomposition at a manageable level, we deliberately fix xmax = 50. It 
has been verified that the same simulations with xmax — 60 generate identical results. 



tively, we can set a maximum value for \. This allows us to obtain results for extended 
time intervals in the range where the system behaviour remains consistent. Figure B~5l 
presents several dynamical quantities obtained from simulating a boson chain prepared 
in a Mott insulator state with one boson per site. The apparent lack of dynamics in the 
graph of the expectation value of the number of bosons can be explained by noticing that 
in a regularly filled chain with one particle in every site the boson hopping generates 
equitable mass exchange, specially around the central part. Therefore, the expectation 
value remains relatively stable through very long times. However, the system develops 
an underlying dynamics, as can be inferred from the graph of fluctuations, given by the 
formula, 
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fi = (^) 2 , (4.14) 

where hi = aja.;. As can be seen, fluctuations saturate faster than the von Neumann 
entropy between complementary blocks. This suggests that the incidence of high order 
effects in the chain is very important, not least in the short run. It also can be seen that 
the system dynamics close to the chain ends is notoriously different from the flat patterns 
displayed over the central part. These border effects are characteristic of chains with 
open boundary conditions |78j . In the long run we expect that these border fluctuations 
propagate all across the chain and the system acquires a wavy profile. Nevertheless, 
long time simulations are likely to contain flawed data as a result of the truncation of 
MPS. 

In this section we have given a short review of the BH model, emphasising in the aspects 
that are more relevant to our investigation. We have described the use of imaginary- 
time TEBD in boson chains and its relation with the phases associated with the model. 
Similarly, we have discussed the main characteristics of real-time numerics and the way 
our simulations arc affected by them. In chains with strong correlations, which are 
the focus of subsequent sections, the presence of entanglement considerably hinders the 
simulation efficiency. As a result, we can estimate the ground state as well as the 
dynamics of highly entangled chains in systems of moderate size. Further analysis and 
complementary results are to be presented in subsequent sections. 
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Chapter 5 

Entanglement and its relation 
with physical processes in 
engineered Bose-Hubbard 
chains 

As the increasing development of cooling techniques in optical lattices has led to the 
experimental implementation of engineered quantum systems, it has become important 
to explore how the quantum nature of such configurations can be used in quantum com- 
putation processing and other interesting fields of research. Consequently, given the 
connection between entanglement and quantum computation, we have concentrated on 
studying the entanglement in systems that are now the object of intensive experimen- 
tal exploration. In doing so we have focused on one kind of entanglement that, to our 
understanding, is most useful for developing quantum computation tasks, this is, the en- 
tanglement between distant places of the chain. Likewise, given the volume of technical 
complexities associated with a multi-disciplinary investigation, we have narrowed this 
research to the study of entanglement between the ends of the chain and the entangle- 
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Figure 5.1: Potential realization of BH model with spatially dependent repulsion in an 
optical lattice. 

ment register provided by the coefficients of the canonical decomposition. In addition 
to the characterization of entanglement, in this section we look at the relation of such 
entanglement with the physical mechanisms associated with the model under scrutiny. 



5.1 End-to-end entanglement and generalities 

In a boson chain with no chemical potential or constant chemical potential on each site, 
the Hamiltonian is given by, 

N u k JV_1 

H = -Y h l h k{dla k - 1) - Moi +1 a k + 4^+0- (5- 1 ) 

fc=l k=l 

The reason for not including additional terms is because this slightly simplified form 
allows us to focus exclusively on repulsion and hopping, so that the effects of chemical 
potential can be considered on the light of the results obtained in this simpler scenario. 
From our initial simulations we observe that ground state entanglement between the 
ends is highly sensitive to the particle distribution profile along the chain. This is 
due to the effect of tracing out intermediate sites in order to get the reduced density 
matrix of the terminals. When population is accumulated in the middle of the chain 
we must trace out states with a considerable number of particles in order to obtain the 
mixed state of two distant sites at the borders. Lower population at the ends implies 
a lower number of possibilities for distinct particle number states available at the ends 
and consequently a lower ability to correlate particle numbers between these sites and 
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Figure 5.2: Left. Particle distribution profile and fluctuations for the ground state 

in PTH chains. ( — ) No repulsion at all. ( ) J72-7 = 100. Right. Ground state 

entanglement between the ends of a chain of eight places and eight bosons against 
repulsion. PTH and repulsionless ends are assumed . Fitting y = 0.91 — 0.34e~ 035x . 



entangle them. Thus for the repulsionless case and PTH, where population accumulates 
primarily in the middle of the chain (which is shown in figure [521, the entanglement is 
vanishingly small. The choice of PTH, as opposed to a uniform Jk, thus clearly illustrates 
the negative effect of population concentration in the middle of the chain. To get a 
higher entanglement, we then proceed to deliberately frustrate boson accumulation in 
the central part of the chain by evenly increasing the repulsion in every site but the ends 
which in turn are left repulsionless (figure l5~Tj) . This approach allows us to look closely 
at the physical mechanisms underlying the interplay between hopping and repulsion and 
its effects on quantum correlations at long distances. In figure I5T21 we can see how the 
boson population in intermediate positions concentrate on the lowest levels of occupation 
creating a fine arrangement of superfluid particles that connects the ends, which in turn 
become highly populated. Figure 15.21 also depicts the behaviour of Log-negativity as a 
function of repulsion. It can be seen that entanglement grows monotonically, undergoing 
saturation for large repulsion, suggesting that maximum entanglement is achieved when 
particles in intermediate places behave as hardcore bosons. In this limit, bosons find it 
difficult to stay at the intermediate positions as the hopping of particles from neighbour 
sites would give rise to multi-occupied states with strong repulsion. Instead, bosons 
organize in superposition states that on average contain only half a boson per site for 
the intermediate sites, creating a globally delocalized profile that favours long-range 
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Figure 5.3: Numerical cost of ground state simulations of Hamiltonian (|5.1[) for N ~ 
M = 20. These computations correspond to some of the points in figures 15.51 and 15.61 
Left. CH. Right. PTH. 



quantum correlations. Similarly, analogous plots corresponding to expectation values 
of number of particles also display saturation for high repulsion, establishing a clear 
signature of rescaling behaviour, that is, bosons reorganize under the influence of larger 
repulsion constants in such a way that the phenomenology of the system is altered, 
roughly, just by a proportionality factor. Similar conclusions can be drawn for chains 
with constant coefficients Jk and repulsionless ends, although PTH chains arc more 
efficient at enhancing entanglement. For instance, in a chain of ten sites, ten bosons, 
PTH and high intermediate repulsion 1/2-9 = 100, Log-negativity was found to be 0.83 
against 0.60 of the same chain with uniform unity hopping. In the study ahead, we 
further explore this remarkable behaviour. 



5.2 Entanglement in the ground state 

The analysis in the previous section suggest that the relative population at the ends is 
in some manner linked to EEE. Indeed, after some test computations we found that as 
a consequence of the universality characteristics of the problem it is appropriate to plot 
our results using the fraction of particles on the ends (dimensionless), 

_ (qjai) + (ajyajy) 

which allows a comparative perspective of the system behaviour for different scales. £ 
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Figure 5.4: Entanglement between the terminals and the rest of the system in a boson 
chain governed by Hamiltonian (|5.1[) and with no repulsion on the ends. Left. CH. 
Right. PTH. As usual, the number of boson in each chain is equal to the indicated 
number of sites. 



varies according to the intensity of the repulsion in intermediate sites. As can be seen 
in figure [5"T2T lcft. when repulsion is small few particles remain on the ends and £ ~ 0. 
From the same graph it follows that when repulsion is strong the amount of particles on 
the ends is maximum and, 



, 1 1 
C =2 + M- 



(5.3) 



Notice that ( keeps a close relation with the repulsion constant Uk- This variable ( 
enables us to depict our results for different chain sizes in a comparative manner. Before 
presenting our results and the corresponding analysis, we want to briefly comment on 
the computational cost associated with the simulations shown on this section. In order 
to make our analysis as robust as possible, we performed computations in chains as large 
as our resources allowed them to be. In a strongly entangled chain, nevertheless, TEBD 
is not very efficient, and the biggest system we worked on corresponds to 20 bosons in 20 
sites. These computations, however, are considerably demanding (using normal desktop 
computers, finding the ground state can take up to two weeks). As a comparison, note 
that the maximum x m the simulations presented in figure l4~4l for systems of 100 bosons 
in 100 places is always below 60 while the same cost-measuring parameter goes well 
above 140 in figure IB~3l which corresponds to the results to be shown in this section. 
In our physical model, independently of the particularities of the hopping profile, the 



87 



ground state of repulsionless chains is eminently superfluid, although with most of the 
tunnelling taking place in the central part, leaving the terminal positions nearly empty 
(figure 15721 - left). In this case entanglement is at the same time strong and localized. 
Indeed, in these circumstances entanglement between half-chain-blocks is strongly en- 
hanced. Similarly, particles are forced to tunnel through longer distances and correla- 
tions develop at longer scales as repulsion in intermediate positions climbs up. Addi- 
tionally, bosons increasingly migrate towards the terminal positions since the repulsion 
prevent particles from staying in intermediate sites. Outstandingly, on account of the 
repulsionless ends, boson fluidity stands no matter how strong the repulsion in inter- 
mediate sites is. As repulsion slowly grows away from zero, entanglement between the 
chain halves goes down while entanglement between the ends and the rest of the sys- 
tem goes up, an indication that the originally localized correlations are being spread all 
over the chain following the particle distribution profile. When repulsion is increased 
even more, the terminals start getting macroscopically occupied so that the expecta- 
tion value of the number of bosons and the fluctuations in intermediate sites both go 
down asymptotically towards J, which indicates that bosons get highly delocalized. In 
this case correlations among places near the ends and in opposite sides of the chain are 
strongly enhanced, in contrast to the correlations in the centre of the chain. This can 
be seen in the graphs of figure [53] where, interestingly, von Neumann entropy between 
both ends and the rest of the system slightly comes down after the original redistribu- 
tion of entanglement mentioned above occurs. This effect is more evident for the CH 
case. The decrease in the entanglement between the ends and the rest of the system can 
be interpreted as a reduction of the correlations in intermediate sites which takes away 
local combination of states with substantial entangling potential. On the other hand, 
this reduction of the correlations between both ends and its complement is convenient 
for the emergence of EEE, since the reduced state of the terminals becomes closer to a 
pure state. This fact alone, however, is not sufficient for the arising of EEE. Even if the 
state of the terminals is completely pure, it is possible that the quantum state of such 
reduced system is separable. In this context it is interesting to explore the mechanisms 
that ultimately define the onset of long-range entanglement. 



ScS 




The issues discussed so far are independent of either hopping profile we feed into Hamil- 
tonian (|5.1[) . We now want to focus on those characteristics associated exclusively with 
the values given to the hopping coefficients. Figure 15.51 summarizes the results we ob- 
tained from computing the logarithmic negativity in both CH and PTH chains. Notice 
that the curves in figure I5~51 are analogous to the results in figure I5~21 which are depicted 
in a slightly different way. 

As a matter of fact, EEE has a more interesting behaviour in PTH chains than in CH 
chains. The most notorious advantage of PTH is the persistence of entanglement, even 
when the number of sites and bosons augment. In contrast, the amount of entanglement 
in CH chains does not stand against increasing chain size. Similarly, when the entan- 
glement values are seen at any given fixed £, it is clear that PTH chains display more 
EEE than their CH equivalent. This is straightforward since in figure |5~51 PTH entangle- 
ment arises well before CH entanglement. Given the increasing cost of the simulations, 
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specially for the PTH case, not all the possible values of £ are considered, but further 
information can be extrapolated from the data at hand. Indeed, PTH curves fit nicely 
into a single line with slope 1.75. As expected, the cost of the simulation is proportional 
to the amount of entanglement, but we can get a reliable scaling characterization from 
the information contained in figure I5~51 Correspondingly, numerical fitting yields 

EEE^ H 5 ~ ;o 5 (7V 01 ),and 

EEE^ 5 ~ log(N-°- 2 ). (5.4) 

Indeed, the logarithmic behaviour of EEE has been inherited from the measure we use 
to quantify EEE, but the fitting formally establishes that the PTH profile provides a 
platform that is highly convenient to entanglement. 

The fact that entanglement is better for some configuration of parameters than for others 
makes us wonder what is so special about a hopping profile that delivers substantial en- 
tanglement between the terminals. Curiously, the particle distribution profile in general 
looks like the one in figure l5.2U eft for chains with very high repulsion in intermediate 
sites, no matter which hopping profile, PTH or CH, we choose. This means that the 
mechanisms giving rise to entanglement must be determined by high order terms of the 
imaginary-time-evolution operator |16j . Such high order terms induce essentially two 
kinds of physical processes, namely, 

• Exchange of more than one particle among places of the chain. 

• Exchange of particles among distant places of the chain. 

These mechanisms can take place simultaneously, and as a result, correlations develop 
at a macroscopic scale. When the conditions are optimal, as for instance with PTH, 
correlations are massively enhanced and this leads to long range entanglement. This 
interpretation coincides with further results, especially with the logarithmic negativity 
graphs shown in figure 15.61 As can be seen, entanglement between the halves of the 
chain decreases monotonically as the repulsion in intermediate sites is turned up. Such 
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of bosons in equal to the indicated number of sites. Top. CH. Bottom. PTH. 
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behaviour is in certain way not surprising since as repulsion becomes more intense, 
more particles are forced away from the centre and into the terminals, which reduces 
the amount of "quantum channels" that can be used by every half-chain block to speak 
directly to each other. Nevertheless, the behaviour of logarithmic negativity in PTH 
chains is quite different. In such a case the graphs display a clear minimum instead of 
a constantly falling curve. Consequently, in spite of the reduction in communication 
resources, there is a critical repulsion value for which entanglement stops falling and 
slowly turns its way up. The bigger N, the sharper the effect, which suggests this 
would lead to a second order phase transition in the thermodynamic limit. This sudden 
reinforcement of entanglement can be understood as coming from an enhancement in the 
particle exchange between the chain terminals. Certainly, when the dclocalization length 
of bosons is equal to the distance between the chain ends, particles held on the terminals 
can be exchanged directly between the ends. This induces a transition of the system 
state, since those particles originally squeezed against the chain borders are suddenly 
unleashed and, instead of being highly localized around the system boundaries, their 
wave functions now envelop the chain to its whole extent, keeping on average most of 
the wave function weight on the terminals. This interesting effect can also be intuitively 
identified by assuming a particle-like behaviour of bosons. Indeed, we can think that 
as repulsion in the middle is continuously turned on, bosons are obliged to hop through 
longer distances and thereby become more and more delocalized. The enhancement 
of entanglement is therefore determined by the tunnelling scope rather than by plain 
particle accumulation, although the latter is a necessary condition for EEE emergence. 
In CH chains, the form of von Neumann entropy between both halves in consistent with 
a regime in which entanglement is being continuously redistributed across the system 
as a result of increasing tunnelling scope, but the fact that there is no turning point in 
figure [5751 -toD indicates that such hopping never takes place across the complete length 
of the chain. However, such hopping is enough to induce EEE at finite N, but not 
in the thermodynamic limit since EEE dies down against increasing chain size. It is 
as if such finite N entanglement results from the overlapping of neighbouring wave 
functions. The transition to a strongly entangled state takes place when such local wave 
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functions flatten and successively combine to create a strongly global description of the 
state. Once the previously mentioned long scale tunnelling takes over, increasing the 
repulsion in intermediate sites reinforces end-to-end exchange and induces an increase in 
the amount of entanglement shared by the terminals, as can be seen in figurc HTBl -bottom. 
This can also be understood in terms of the so called healing length, which defines the 
distance that takes the condensate to form its bulk. We can say that the healing length 
determines the space over which the superfluidity of the system is suppressed. In terms 
of physical constants, this length scales as Lhealing ~ ~Ju' ^ n 4 -^ e ! the healing length 
is around O.lnm. The phenomenon that we report in this section can be seen as a 
suppression of the healing length, in such a way that superfluid features are present all 
over the chain, and then these features are reinforced by the PTH hopping profile. 



5.3 Entanglement and dynamics 




Figure 5.7: Logarithmic negativity for N = M = 8 and CH. The system is prepared 
in a Mott insulator state with one boson in each site. As can be seen, strong repulsion 
constants in intermediate places are necessary to induce substantial EEE. 
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Here we study how entanglement arises dynamically in boson chains with superfluid 
characteristics. This is the most relevant situation since any form of entanglement is 
null in the Mott insulator state. In the discussion that follows we assume that our 
chains arc initially prepared in a Mott insulator state and the parameters are suddenly 
changed to generate the dynamics. The new set of parameters correspond to a given 
hopping profile, CH or PTH, and a given repulsion profile which always maintains the 
ends repulsionless. This is the most natural scenario one can expect in a BH chain 
where one first prepares the Mott state with one boson per site and then turns all the 
interactions globally, the same as in [48 , 79 , 80] , but opposite to the approach in [811 [78] 
where the quench goes from the superfluid to the Mott insulator. As a generic example 
we present in figure 15.71 the evolution of the logarithmic negativity for different values 
of intermediate repulsion. In a chain with all the repulsion constants set to zero, EEE 
shows very little development in the course of evolution as can be seen in the indicated 
figure. Dynamics in repulsionless chains is, in fact, known to lead to a progressive 
thcrmalization of reduced density matrices of each site [791 180j . Two site entanglement 
is thus very low in this repulsionless scenario. High repulsive evolution, on the other 
hand, shows a less regular profile, product of interference among many unsynchronized 
phases accumulated locally. From figure 15.71 we infer that as in the case of ground state 
entanglement, dynamic entanglement does not build up unless a significant accumulation 
of particles assists the enhancement of fluctuations at the ends. As can be observed in 
the evolution patterns of chains with high repulsion in intermediate sites, the amount 
of entanglement contained in the system after reasonable times can be significant, yet 
smaller than what can be obtained from the ground state in optimal configurations 
(compare with figure [5~2l right) . Here, we assume uniform hopping but our conclusions 
are valid also for PTH. 

In order to identify relevant differences between the entanglement behaviour in chains 
with different hopping profiles, we found it necessary to perform simulations in chains of 
up to 20 places with 20 bosons. As we already know that strong repulsion is necessary 
to induce EEE, we just set high repulsion constants on every site except the terminals 
and simulate the evolution using the TEBD algorithm. The resulting dynamics can be 
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Expectation Value of the Number of Particles 



Expectation Value of the Number of Particles 




Figure 5.8: Particle distribution in bosonic chains for N = M = 20. The system is 
prepared in a Mott insulator state with one boson per site. Strong repulsion constants 
in intermediate sites forces the bosons into the ends. Left. CH and C/2-19 = 100-Er. 
Right. PTH and U 2 -w = 200E R . 

seen in figure 15. 81 where the expectation value of the number of particles is plotted as a 
function of time. The ends get macroscopically occupied and tunnelling in intermediate 
places intensifies. EEE arises some time after the ends have been occupied, as can be 
appreciated in the graphs of figure IS~9"1 In general, the bigger the chain the longer it takes 
for EEE to emerge. For PTH, a natural time scale is determined by the transmission 
period, T = tt, on account of our particular choice of A = 2. Entanglement shows 
up more or less at half of the period, t « when some particles have travelled from 
one end of the chain to the other. It could be, however, that EEE in longer chains 
behave differently, since repulsion in intermediate places hinders the boson mobility. 
This negative effect is likely to be more pronounced when transport takes place at long 
distances. For longer times entanglement might grow above the values reported in figure 
15.81 independently of the hopping profile, but this effect is hard to follow as long time 
simulations require additional computational resources |77l I76j . 



5.4 Entanglement and perturbative dynamics 

We have seen how boson chains can be made to contain substantial entanglement in 
the ground state as well as in the course of dynamics. Ideally, if both approaches 
are combined we then could multiply the amount of entanglement in the chain, which 
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Figure 5.9: The system is prepared in a Mott insulator state with one boson in each 
site. These graphs correspond to equivalent systems shown in figure [HTSl Left. CH and 
U 2 -ig = lOOEj?. Right. PTH and U 2 -ig = 200E R . 



would result in a highly beneficial use of the several quantum degrees of freedom that 
come from the multiple possibilities of arranging a given number of particles on the 
system. Hence, we now assume that a chain initially prepared in an entangled ground 
state evolves by the action of a perturbation added to the original Hamiltonian. The 
reason for choosing a perturbation rather than an abrupt change of parameters in the 
Hamiltonian is because a perturbation is expected not to modify the particle distribution 
profile substantially, so that the terminals of the chain remain macroscopically populated 
which has been already identified as a necessary condition for entanglement generation. 
In formal terms, the state evolution is given by, 



m))=e- itH ^ g ), (5.5) 

where. 

N 

H 1= H + eJ2hj^ (5-6) 

in such a way that \ip g ) is the ground state of Hamiltonian Hq 7 from equation (15.11) . and e 
is a small real number which determines the intensity of the perturbation, hj represents 
a local operator acting on site j. The explicit form of operators hj is to be worked out 
according to our convenience. For consistency, the first condition to be satisfied is, 
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Figure 5.10: Sketch of bosons spread across a chain of repulsionless ends. 



[H ,J2hj]^0, (5-7) 

otherwise dynamics would be trivial. In order to simplify our analysis we assume that the 
local operators can be written as functions of the corresponding local number operator, 



hj = hj(fij) (b.- 

where, 



hj = OjCLj, (5-9) 

so that we can now focus on the functional forms hj (x) rather than in quantum opera- 
tors. We can establish a non-vanishing commutator between the perturbation and the 
Hamiltonian by defining, 



hj{x) = kx, (5.10) 

where k is the integer distance between site j and the closest end, as indicated in 
figure 15.101 This form holds only for intermediate places. Similarly, to determine the 
form of the operators on the terminals, /ii,./v, the primary factor to take into account 
is the optimization of boson-transfer from intermediate sites towards the ends. If the 
expectation value of the number of particles remains stable during time evolution we 
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can carry out a classical analysis in the basis of the number of bosons. Let us think that 
in the configuration shown in figure 15.101 which sketches one of the particle distribution 
shown in figure [5~2l left, we take half a boson from the centre of the chain and put it in 
the closest terminal. We can assume that this repositioning docs not cause much change 
in the average value of Hq, since there is no contribution from the repulsion terms for 
bosons in intermediate positions, because repulsion is zero for one or zero bosons, or for 
bosons on the chain ends, because repulsion is off on the terminals. As a consequence, 
we would expect that the change in the average values of the system could be calculated 
focusing only on the averages of the perturbative potential. Additionally, because we 
want particles to go from the middle of the chain to the terminals and vice versa, we 
need that the energy cost of extracting particles from intermediate positions balances 
the energy cost of putting the same particles on the chain ends. For example, starting 
from the state depicted in figure 15.101 if all the particles go consecutively from the 
middle of the chain to the ends, the first process to take place is that bosons in positions 
with k = 3 disappear from their respective positions and pop up on the terminals. The 
corresponding energy balance can be written as, 



In this expression no represents the expectation value of the number of particles on one 
of the terminals, which can be simply taken as no = 2.5 for this specific case. Once 
this process is over, there are fewer particles in intermediate sites and more on the ends. 
Hence, taking first bosons with k = 2 and putting them on the ends and then doing the 
same for bosons with k = 1 we obtain, 




(5.11) 




(5.12) 



In a chain of size N these equations lead to one single expression, namely, 
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/i ( no+ K^ _fc ) + 0" /i ( no+ K^" fe )) = ¥- 



(5.13) 



After some preliminary tests we find that the optimal choice for the unknown function 
is, 

h(x) = c 2 x 2 + c\x, (5-14) 

Constants c 2 and c\ are used to balance equation (|5.13[) . Combining these equations 
and after some simplification we get, 



c 2 (no + Hf "*)) + ? + f = f 

^c 2 7i + ^ + f + f-^ = f (5.15) 



In order to work out a set of constants valid for every k we must define, 



C2 = -e, 

Cl = (^P + 2n °) (5 - 16) 

Moreover, for a perturbation written as above, the mean number of bosons does not 
fluctuate much since the dynamics is still governed by a Hamiltonian with strong re- 
pulsion coefficients in intermediate sites and bosons are forced to stay on the ends. In 
this situation one would expect that the evolution is dominated by high order dynam- 
ics such as the exchange of particles between the terminals and the rest of the system. 
Evidence of this can be seen in the graph of fluctuations shown in figure 15. Ill -left. As 
can be expected from our particular choice of function h{x), the transit of bosons be- 
tween the ends and the rest of the chain is enhanced and therefore the fluctuations on 
the terminals are developed according to the strength of the perturbation. Certainly, 
such enhancement of fluctuations leads to an improvement of communication between 
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system responds according to the intensity of the perturbation. Equivalent graphs cor- 
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Fluctuations in the number of bosons allocated in one end of a chain. Right. Dynamical 
entanglement. 



distant places which results in entanglement. Figure [5 . 1 U right shows Log-negativity as 
a function of time for different perturbation intensities. Roundoff errors accumulated 
in the ground state during the imaginary time evolution cause slow oscillations in the 
dynamics at e = 0, but these are tiny compared with the more complex evolution seen 
at finite e. In the initial stages of evolution the dynamics is characterized by a sharp 
increase of quantum fluctuations on the ends of the chain accompanied by little change 
in the expectation value of the number of bosons. Hence, entanglement development is 
enhanced while avoiding massive migration of bosons towards the centre of the chain. 
This effect is maximal at around e = 0.1 and becomes less important for higher values 
of e, where perturbative terms induce a less predictable dynamics. Significantly, entan- 
glement generation is improved not by adding interaction terms to the Hamiltonian but 
by evolving an already entangled state using a locally-tuned perturbative potential. 
In order to see if this effect lingers on when the size of the chain is increased, we 
applied the same ideas to chains of N = M = 14. Interestingly, EEE behaves in a very 
similar way than in the previous case, with almost the same values at the same times, as 
indicated in figure [5T41 In this figure we also include CH results, which make it clear that 
the enhancement effect is characteristic of PTH chains. This is a curious matter, since 
the perturbative analysis does not include any reference to any specific hopping profile, 
it just balances the energy costs coming from different terms in an optimal trajectory. 
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Figure 5.12: Logarithmic negativity in chains of N = M = 14 for both CH and PTH. 
PTH lines are thicker than their CH equivalents. EEE in the ground state can be 
improved as a result of dynamics in PTH chains. \ = 200 for PTH and \ = 100 for CH. 



In principle, one can think that the lack of entanglement development may be due to 
either, 



• Lack of boson mobility. 



• Dynamics with deficient entangling capability. 



In general, entanglement seems not to team up with CH, as opposite to PTH. What- 
ever the specific reason why this compatibility between entanglement and some specific 
profiles occurs, our results clearly suggest that particle transport plays a crucial role in 
the development of correlations in the quantum state. 
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Figure 5.13: Entanglement detection 

5.5 Entanglement detection 

In this short section, we want to explore how this entanglement could be measured in 
the laboratory. Once bosons are cooled down in the ground state or after evolution 
has taken place and the interactions are turned off, the detection protocol studied in 
[521 153"] could be used. The idea behind most entanglement- verification procedures is to 
uncover the non-local behaviour of the entangled state. This non-locality can manifest 
itself in a spatial basis, such as it happens, for instance, while probing the non-locality 
of a single particle |56j . but it can also arise in the basis of the number of particles. 
The protocol proposed in [52] exploits the idea presented in [55], that there cannot be 
a pure separable state with fixed particle number and full interference among spatial 
modes. Let us suppose that after the intermediate degrees of freedom in the chain has 
been reduced, we are left with a density matrix p, which represents our knowledge of the 
resulting two-side mixed state, or equivalently, the state of the chain ends. Following 
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the protocol of reference [52], bosons on the terminals are sent through a 50:50 splitter 
and then the number of particles on the outputs is measured, for example following the 
technique proposed in [T5] . 

The effect of the 50:50 splitter is to combine the modes of the original system, d\ and 
ajy, according to the following transformation, 



4 = ^(4-4), (5.17) 

where a\ and a\ are the modes of the system after the 50:50 splitter. Such splitter can 
be considered as the physical place where interference between the bosonic bulks takes 
place. In this way, as a\ and a c are the operators describing the particle occupation 
on one of the outputs, EEE can then be detected from the experimentally measurable 
variable, 



N + 

£ab = tr(ala c p) - — - = tr(a[a N p), (5.18) 



where, 



tr{ala c p), (5.19) 

is the number of particles in one output and p is the reduced density matrix of the ends. 
As it is characteristic of a good entanglement measure, 

e AB = 0, (5.20) 

for separable states and, 

e AB > 0, (5.21) 
for entangled states. In this way, entanglement can be detected by just comparing the 



103 



amount of bosons in one output, which can be done with relative efficiency. In order 
to show how this method can be applied to our model, we present in figure 13751 a graph 
of €ab obtained from the ground state of chains with PTH and CH for some of the 
parameters in figure 15.51 Wc conclude that cab correctly determines whether the state 
is entangled or not, but it does not faithfully depicts the actual amount of entanglement 
in the system as it is shown in figure 15.51 This might be due to the fact that eab 
accounts for the entanglement contained in the degrees of freedom corresponding to first 
order physics, that is, those given in terms of averages and mean values of measurable 
observables, while most of the entanglement in the system comes from high order physics, 
as for instance long range particle exchange. In any case it would be of great interest to 
investigate this entanglement measure in the laboratory, since this would definitely tell 
whether long range entanglement can be achieved in bosonic chains. 
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Chapter 6 

Kicked dynamics in systems of 
ultracold atoms 

6.1 Dynamics of a double-kicked bosonic condensate 

So far we have studied diverse aspects of bosonic systems described by fully quantized 
Hamiltonians. Such an approach is necessary when we want to deal with quantities 
such as entanglement, since one needs to treat sections of the system as individual 
entities. For instance, in order to find EEE we need to get the reduced density matrix 
of the chain ends, which describes the statistics of such two-site system. As we have 
shown, this description of the state bears some technical complications. For this reason, 
alternative approaches are of interest. In cases when the object of study depends on 
the system's behaviour as a whole and not on the combination of mutually interacting 
contributions, a simpler approach provides valuable insight. It is in this context that 
we want to address the problem of the kicked condensate [HH E31 [S3] . Here we want to 
study the stability of a condensate of bosons when it is subject to a periodic kicking. 
The way we manoeuvre the problem is the following: we first study the evolution of 
the non-linear Schrodingcr equation and find the condensate wave function. Then, we 
treat this solution as the ground state of the system, and analyse how the condensate 



105 



population migrates towards excited levels. It is worth mentioning that this latter 
part is carried out following the procedure introduced by Castin and Dum in reference 
[85) . This procedure is in many aspects similar to the Bogoliubov method to study 
interacting Bose gases. However, in the latter case the state of the system is described 
by a coherent state and therefore the number of particles is not well defined. The 
Castin and Dum approach overcomes this issue by expanding the atomic field operator 
in a sum of condensed and non-condensed parts and assuming that the non-condensed 
contribution can be considered as a perturbation with appropriate scaling behaviour with 
respect to the condensed part. In this way, the total number of particles is always fixed 
and the proportion of particles in the condensate can be calculated as the expectation 
value of the condensate field operator. Suppose we have a bosonic system which can 
be described by means of the Gross-Pitaevskii equation, which corresponds to a mean 
field approximation of the Hamiltonian (|4.1[) . plus a periodic kicking. The particles are 
confined in a ring-shaped trap of radius R. We assume that the lateral dimension r of 
the trap is much smaller that its circumference, and thus the system can be effectively 
considered as one-dimensional. In this model, the one-particle wave-function, "0, is given 
by the following second order differential equation, 



^ = + SMV + kcos6 (f^S(t- n/3) - S(t - n0 - e) W. (6.1) 

\n=0 J 



This model can be physically associated with a Bose Einstein Condensate (BEC) that 
is being double-kicked with periodicity j3. Two consecutive kicks have opposite sign so 
that they nearly cancel each other and depletion is contained. The time interval e is 
a small number and can be considered as a perturbative parameter. The intensity of 
the external potential is k and the interaction strength is g, which is described by a 
mean-field term with scaled strength [531 IM] , 

8Ma s R 

9 = r2 , (6-2) 

where a s is the s-wave scattering length, and M is the total number of bosons. The 
length is measured in units of R and energy in units of, 
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with m the atomic mass. In addition, we assume periodic boundary conditions ij>(0) = 
ip(2n). It is worth mentioning that in writing equation (|6.1[) we have discontinued the 
notation of previous sections in order to make our formulae compatible with standard no- 
tation in the field. Equation (|6.1[) is non-linear, hence conventional quantum-mechanics 
methods cannot be applied. In simple terms, our problem consists in getting an analyt- 
ical solution of the GP equation for the initial condition, 

M8) = 4=- (6-4) 

Let us point out that because this initial wave function is symmetric with respect to 9. 
the wave function must reflect the same symmetry for all times, since such symmetry is 
explicit in the GP equation as well. Equation (|6.ip along with the initial condition (|6.4[) 
can be solved straightforwardly for k = or e = 0. The case g = is a double- kicked 
quantum rotor [861 1871 188] . The time evolution of the initial wave function is dictated 
by the Floquet operator, 

F(f3,e) = f 1 {l3)e %kcoa[e) f 2 {e)e- lkcos[e \ (6.5) 

where T\ ((3) and T2 (e) are the evolution operators for times j3 and e respectively. As T2 (e) 
induces evolution for a small interval of time, the wave function can be approximated 
by its average value during this short period so that we can write, 

This expression is exact to first order in e. As can be seen, the right-hand side of the 
equation is reminiscent of a unitary transformation. It is not difficult to see that such a 
transformation induces the following changes in the quantum operators, 



o9 ou 
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e ikcos W\iP fixed (6)\ 2 e- ikc ° s W = \^ f . lxed (0)\ 2 . 
Inserting these results in equation (|6.5[) we obtain, 



(6.8) 



F(8,e) « fi(/3)e ie (^^ +lfesin ^ 2 " 9|,A/ " e<i|2 ). (6.9) 

When written in this way, the Floquet operator allows us to identify the effective per- 
turbation due to two consecutive kicks. It is worth recalling that efc must be small so 
that the net effect of two consecutive kicks can be considered as a perturbation and the 
wave function evolves slowly as seen stroboscopically after pairs of kicks. Similarly, it 
should be noticed that equation (|6.9p is exact for g = 0. We would like to find a way of 
introducing a perturbative potential in equation (|6.ip from the information contained 
in equation (|6.9[) . Curiously, should the expression on the exponential be exclusively 
^-dependent, it would be easy to carry out such a procedure. We then could simply 
introduce a perturbative potential like, 

oo 

SV(9,t) =f(6)J2*(t-M- (6-10) 

However, the operator on the exponential in equation (|6.9[) depends explicitly on ^ and 
therefore dynamical effects take part in the perturbation process as well. To see this 
more explicitly, let us put the term on the exponential in equation (|6.9[) in the following 
form (up to an i factor), 

{ I W ~ e9 ^f ixed \ 2 } + % \ k sin ^ + l \ k cos B ~ | fc2 sin2 e - (6 - n) 

The term in brackets represents free evolution while the other terms account for the 
perturbation. We would like to disentangle these contributions in order to simplify 
the dynamics. As the commutator between terms is second order in e, splitting the 
exponential preserves first order accuracy. Similarly, it can be shown that as long as 
we stay in the small perturbation regime, the term that goes with the first derivative 
actually depends on e 2 , hence we drop it. After these simplification the Floquet operator 
can be written as, 
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F(fi,e) « fi(/3)e ie (^^" fflv ' /la!e£!|2 )e i ( i 4 fccos ^"4' i:2sin2 ^. (6.12) 
Furthermore, putting together terms representing time evolution and obtain, 

f 1 08)e fc ('^ r_9l * /i " d|a ) = T(/3 + e). (6.13) 

This step amounts to sticking the short time slice between kicks to the longer time 
interval between pairs of kicks. This sticking is mathematically correct and does not 
introduce additional errors to our analysis except by the fact that ip is considered time- 
independent in the short time interval. 

Introducing expression (|6.13[) in (|6.12l) we hnd that the Floquet operator can be written 
as the kicking term followed by free dynamics. Moreover, from the Floquet operator we 
can go back to equation (|6.1[) and rewrite it as, 



i^=-^+^|V+(-4 C0S W + £ r sin2 w) (£*(i-n(r))W, (6-14) 

where T = [3+e. Now the label n in the sum runs over pairs of kicks rather than on single 
kicks. Correspondingly, the original problem has been re-formulated in terms of better- 
understood single-kick dynamics [HH HOI GED ISZ] • The fact that the term proportional 
to cos(#) is not hermitian does not lead to any inconsistency because the whole term 
accompanying the delta function does not show up in the energy formula (equation 
(|6.36[) ahead). This complex contribution in the new GP equation is related to kinetic 
effects coming from differentiating equation (|6.9p . 

In what follows, we apply the Castin and Dum linearisation-approach presented in ref- 
erences [HU |55] . Briefly, the formalism consists in linearising the time-dependent GP 
equation for a small time-dependent perturbation. Such approach leads to a set of two- 
variable differential equations driven by inhomogeneous terms. The condensate evolution 
is then determined by the form of the small time-dependent perturbation or source. In 
order to review the Castin and Dum's method [52] , let us write the solution of equation 
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(|6.14[) as a sum of two terms, 



il>(9,t) = + 81>(0,t), (6.15) 

V Z7T 



so that replacing (|6.15D in (16. 14[) we get, 



where, 



.dSi> 1 dH^p 2gSj> gSr SU 

= -vie + — + ^r + v^> (6 - 16) 



SU{9,t)= ( -iy cos 8+ sin 2 o\ ( 6{t - n(T)) \ . (6.17) 



\n=0 

|2 



In addition, we have neglected a term proportional to |<5'0| 2 m equation (|6.16[) . which 
means we have removed the non-linearity. As the wave function is complex, the real 
and imaginary part can be treated as independent real functions. Nevertheless, it is 
more elegant to keep the original function Si]) and introduce its complex conjugate as 
the other unknown. It can be shown that in doing so equation (|6. 16|) leads to, 




l 

2 W 1 



<J 

2jv 



a 

2-7T 




(6.18) 



In our case, the source term is given by, 



S(0,t) =5U(d 1 t)2jj = (-^cos# + ^-sin 2 0) lJ26(jb-n(T))\-±=. (6.19) 



2 2 J \^ v ' " V2n 

' \n— / v 

This is essentially a linear equation and as such it can be solved using standard methods. 
Some care must be taken, for instance, in order to make sure that function Sip has a 
physical meaning so that it remains orthogonal to the homogeneous solution ipo = 
Issues like this are extensively discussed in reference [92.. In what follows, we limit 
ourselves to using the results of this reference to get a solution to our problem. In 
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synthesis, the wave function is given by, 



W,t)= (6.20) 
1 (1 + UMty 9 + Vxbl{t)e- i0 + U-xb-x (t)e~ ie + V-^tye" 



V2ir 

+U 2 b 2 (t)e 2W + V 2 b 2 (t)e~ 2W + U- 2 b_ 2 {t)e~™ + V- 2 b*_ 2 (t)e M8 ) . 
The unknown coefficients can be obtain from time integration, 



Jo 1 



where, 



(6.21) 



Sj (t) = JJdB {Uj S(9,t) + Vj (^L\ S*(0,t)) , (6.22) 

and finally, Uj and Vj are given by [Hlj , 



U 3 +V 3 =Q, 

Uj - Vj = 1 (6.23) 

with, 



and, 




(6.24) 



(6.25) 



In principle, j can be any integer number. Fortunately, only four modes contribute, 
namely, j = —2, —1, 1, 2. Performing the integrations we find, 
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*i(*) = ^(-Ui + Vj !f^6(t-n{T))\ , (6.26) 



,2 / \ 

s 2 (t) = ~ e —(U 2 +V 2 )(J2 S(t - n(T)) , (6.27) 

and s_i = si, s_2 = ^2- Now we must calculate the time-dependent coefficients (|6.21[) , 
which are found to be, 



h(t) = f (-l/k +^i)E^io" 1 e- 8e l(t "" T) , (6.28) 

fc(f) = ^(^2 + K 2 ) E^"* e-"»(*-» T ), (6.29) 

6_i(t) =6i(t), (6.30) 

b_ 2 (t) = b 2 (t). (6.31) 



In equations (|6.28[) and (|6.29D iV(i) is the number of pairs of kicks (N(t) — 1, 2, ...) and 
£1,2 are the eigenenergies associated to each mode. As can be seen these coefficients are 
written in terms of time, nevertheless, it is more convenient to put them in terms of the 
number of kicks in order to compare with numerical results, which lead us to replace t 
by (N — 1)T. After some rearrangement coefficients b\ and b 2 adopt the following form, 



iei(N 



N-l 
n=0 



4Ci sin 



2 



(6.32) 



b2{N) = i e Jfk e -« 2 (JV-i)T y e « 2 T„ = • efc 2 C 2 -i^ {N -i) gin ^jv, (6 33) 
8 fo 8sin%- 2 



On the other hand, the wave function symmetry can be made explicit by factorizing the 
exponentials in equation (|6.21|) . This leads to, 
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V>(0, N) = (l + 2(Uih(N) + VxbKN)) cos9 + 2{U 2 b 2 (N) + V 2 b* 2 (N)) cos 26) . (6.34) 



Furthermore, inserting 61 and b 2 in ij)(8, N) and after some direct simplifications we 
obtain, 



ip(6,N) 

V2tt V 2smwi V Ci 

efc 2 



1 sinwiTV coswi(iV — 1) 5- sinwifiV — 1) cos 6 

°3inu;i V Ci / 

cj 2 N ((% smu 2 (N - 1) + i cosw 2 (iV-l)) cos 26>^ , (6.35) 



■ sine 



4 sin W2 

where w 7 - = . This wave function reads —h= for e = or k = as expected. In order 
to check the accuracy of equation (|6.35p we calculate the state energy by the formula, 



E(N) = J n der(N) (-\^ + \g\^ Wl 2 ) m), (6.36) 

and compare it with computer simulations that use fast Fourier transform to shift the 
state wave function from position to momentum representation and vice versa. The 
graphs in figure 16.11 show in detail the energy oscillations for different values of k and 
g both from numerical simulations and from wave function (|6.35[) . Outstandingly, both 
approaches agree quite well over a wide range of parameters. 



6.2 Stability Analysis 

So far we have considered the condensate as an isolated system. However, in a more 
realistic model particles in the ground state interact with particles in levels of higher 
energy and this causes bosons to escape from the condensate. Both the repulsion as well 
as the constant kicking determine the system stability. The second quantized analogous 
of the GP equation is given by [55], 



113 



Energy 



Energy 




40 60 
Pairs Of Kicks 



LnciL'y 



0.16 

0.1598 | 

* 

0. 1596 | 
111 

0.1594 
0.1592 



llillilSi 



40 60 
Pairs Of Kicks 



LlKTLr)' 




40 60 80 

Pairs Of Kicks 




Figure 6.1: Right graphs show results obtained from equation (|6.35[) and left graphs 
show the same results obtained numerically. Top. k = l,g = 1. Medium position. 
k = 1,9 = 2. Bottom, k = 2,g = 2. Energy is measured in terms of the recoil energy 
given in equation 



H = 



(6.37) 



where, 



W(6,t) 



ek ek 2 

-i — cos(9H sin 9 

2 2 



(6.38) 



^sin 2 c9 £ 5{t-n{T)) 



In a perturbative approach the field operator can be written as, 
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il>(9,t)=il>{9,t)ao+5il}(6,t), 



(6.39) 



where ip(0,t) is the ground state wave function of the condensate, given by equation 
(|6.35|) . &o and Sip{9, t) are the bosonic operators for particles in and out the condensate 
respectively. Hamiltonian (|6.37p must be inserted in the Heisenberg equations of motion 
and then these equations must be analytically integrated to get the Heisenberg dynamical 
operators. Ideally, the number of non-condensed particles (NNP) can be obtained by 
adding up the number of particles in excited modes. However, this can be rarely done 
directly due to the non-linearity of the original Hamiltonian. For this reason we use 
the method of reference [85] , where practical equations in the linear response regime are 
proposed for the non-condensate cloud. From this work we now borrow the following 
identity, 



A(0,i) 
At(0,t) 



Vj(6,t) 



' V*(9,t) 



> - 



(6.40) 



where, 



A = 



(6.41) 



Although A is not the non-condensate field operator itself, it is relevant to the number 
of particles off the condensate since, 



(6.42) 



in the low temperature limit [85) . The number of excitations for each mode in equa- 
tion (|6.40|) is completely determined by the temperature through the well known Bosc- 
Einstein distribution, 



1 



(6.43) 



where /3 = jr-?- At zero temperature (bjbj) = and NNP is dictated only by the 
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components of Uj(9,i) and Vj(9,t) in the space orthogonal to ip(9,t). The evolution of 
Uj(6,t) and Vj(6,t) is dictated by, 




(6.44) 





' w 


^ 




)• 










V 


-W j 






with VF = and -0 = ^(M)- Additionally, NNP can be written as, 

NNP = J2(V j \QQ\V j ), 



(6.45) 



with, 



Q = i- 



(6.46) 



in the limit of very low temperature. In order to carry out our numerical analysis, we 
first integrate the GP equation (|6.1[) using a fast Fourier transform algorithm and get 
the wave function ijj(9,i). Then we insert this ijj(9,i) into equation (|6.44[) and integrate 
for the initial conditions (|6.66[) . Finally, we calculate NNP using equation (|6.45[) . 



6.3 Evolution of the number of non-condensed parti- 
cles 

Equation (|6.45[) is the same as 

NNP = J2{{V j \V j )-\(m)\ 2 }- (6-47) 

3 

Moreover, as ip is normalized to 1, NNP is always positive and dependent on (Vj\Vj). 
Correspondingly, we focus on finding the norm of Vj(9, t) since it provides us with quality 
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information regarding the NNP. For g = 0, Uj(8,t) and Vj(9,t) become decoupled in 
equation (|6.44[) and both functions undergo unitary evolution that leaves the norm 
unchanged. This implies that norm diffusion is exclusively due to the non-linear term 
of equation (|6.44[) . In order to explore how the norm of Vj(6,t) behaves in time, we 
highlight how the non-linear term acts on the state vector, namely, 



where, 



Uj(6, t + A) 

Vj(e,t + A) 



I 



JAG(8,t) 



U 3 (0,t) 
Vj{B,t) 



(6.48) 



G(6,t) = 



2 5 |V(M)P 



-gij}{e,t) 



-2ff|V(^*)| 2 



(6.49) 



As can be seen, operator G is not hermitian, which is by no means inconsistent, since 
further exploration shows that the eigenvalues of G are themselves real. This is an 
important issue to consider when exploring the development of NNP. In fact, according 
to [HI], the norm of the Uj-Vj vector (global norm) is given by, 



[ \Uj) 


H 







1-k rl-K 
2 



{(Uj\, -{Vj\) \ \= de\Uj(e,t)\ 2 - / d9\V 3 (6,t)\ 2 = Constant, (6.50) 



with Constant = 1 for all j. As can be seen, the norm of Vj{9,t) can grow indefinitely 
as long as the difference of the two norms remains constant. Notice that norm growing 
depends heavily on the definition of inner product in the Hilbcrt space. If we had + 
instead of — above the norm of Vj(0,t) and Uj{6,t) would be bounded by Constant. 
Making use of norm conservation properties we now write the state vector as, 



(6.51) 



' Uj(e,t) \ _ I sinh («,-(*)) «,•(*,<) 
v Vj{9,t) j \ cosh( aj (i)K(M) 

ij(6,t) and Vj(6,t) are complex functions normalized to unity while aj(t) is a real 
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function that characterizes the norm of Vj(6, t) and Uj(8,i). In order to gain insight 
into the evolution of the norm of Uj(9,i) and Vj(8,t), or equivalently, the evolution of 
a(t), we explicitly perform the operation indicated in equation (|6.48|) following a similar 
procedure as in reference |93| . In doing so we find that matrix G can be written as (let 
us omit the dependence with 9 and t for simplicity), 



G = 



2<#| 2 sMV* 



with the following eigenvectors (not normalized), 



(6.52) 



(e\E+) = 



(0\E-) 




\+ = V3g\iP(6,t)\ 2 , 
A_ = -V3g\^(d,t)\ 2 . 



(6.53) 



(6.54) 



From equation (|6.48D we can infer, 



V 



t+A 



(6.55) 



Scalars fj,j(9,t) and Vj(9 1 t) can be written in terms of Uj(9,t) and Vj(9,t) making use 
of this equation with A = 0. In doing so we find, 



M(il) = 1 4jM^W| i (6 56) 

Vjm= j^±mm±£^MA. (6 . 57) 

^ ' 2(3 + 2^3) V ' 



At this point it is worth noticing that global-norm conservation is a direct consequence 
of, 
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(E_\E+)=0. (6.58) 

Now, let us introduce the symbol [...|...] to describe an inner product as in equation 
(I0.50|) but with + instead of — , so that we can address the following "global- norm" , 



((U j \U j ) + (V j \V j )) t+A = (6.59) 
([E + \\ H \ 2 \E + ] + \E_Wvtf\E_] + [E + \^e-^+ A \E_} + [£», e 2 * A + A |£ + ]) t , 

where use has been made of equation (|6.55[) . If we expand exponentials to first order in 
A we obtain, 



{{Uj\Ui) + {V 3 \V 3 )) t+A = {{Ui\Uj) + (Vj\Vj)) t + (6.60) 
2*A ([E-\\+v]n j \E+] - [E + \\ +f i*v j \E_]) t , 



hence it follows, 



{JJJjPjl + (V 3 \V 3 )) t+A - m\Uj) + {Vj\Vj)) t 
A 

2i {[E_\\ + v* N \E + ] - [E+\X +f j,*u j \E_]) t , 



which amounts to, 



(6.61) 



-8V3g / de^Imiv*^). (6.62) 
Jo 

So, using (|6.56[) and (|6.57D to write the integral above in terms of Uj and Vj we get, 



2tt 



j t {{UiPi) + tyj\Vj)) = -*g jf dei^Imie^VjU*). (6.63) 
This expression can be further simplified by writing it in terms of aj(t) as defined in 
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equation (|6.5ip . 



d 



(sinh 2 OLj + cosh 2 otj) = — 4<jsinhaj cosh otj / d9Im(ip 2 u*Vj). (6.64) 

Jo 

Finally, using well known identities we arrive to 



= -g J d9Im(ilj 2 u*Vj 



(6.65) 



This equation determines the evolution of the state norm in our problem. Let us point 
out that this equation is valid for arbitrary t since both the dynamical term as well as 
the static potential in equation (|6.44|) do not affect the state norm. Now equation (|6.65D 
can be used to understand some aspects of NNP evolution. Indeed, for t — all the 
variables in the integral are known, namely, 



1 e ij 

Vj = 7^' Uj = 7^> Vj = V^> (6 - 66) 

from which it follows, 



daj 
~dt 



= 0. (6.67) 

t=o 



Consequently, as long as the kicking is not strong enough to drive the condensate and 
the modes out their initial states, the NNP will remain constant (not necessarily zero). 
This behaviour is positively corroborated by the graphs that we show in figure Wl2[ left, 
where the NNP is plotted against the number of kicks. As can be seen from this graph, 
dissipation is quite moderate even after several tens of kicks. 

Figure 16.31 shows the behaviour of NNP for different values of the parameters g and k. 
Initially, NNP does not respond to the kicking, which corresponds to the time that the 
system takes to evolve away from the initial conditions, in which case particle diffusion 
is given by equation ()6.67[) . Once this initial unresponsiveness is over, NNP starts 
to develop in a characteristic manner. In the small interaction regime, NNP-growth 
becomes linear with a slope equal to 2. Empirically, this dissipation is tolerable and 
the condensate can be considered as stable [521 IHU- In contrast, for g = 5 NNP-growth 
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Figure 6.2: Left. j3 = 1.96, e = 0.04. It is seen how NNP is hardly altered for weak 
kicking. Right, fc = 1, <jr = 1. Momentum distribution for t/j and u\ after 160 kicks (80 
pairs of kicks) . After several kicks the momentum distributions still look very similar to 
the initial ones. The condensate and the modes develop opposite phases and the integral 
in equation (|6.65[) vanishes. 



becomes exponential, which unequivocally characterizes the unstable phase. Figure 16. 2X 
right shows the slopes of several curves of NNP as they are all fitted to exponential 
behaviour. Although exponential behaviour may not be the most suitable fitting for 
every curve, it allows us to put our results in a comparative perspective. Figure 16.31 
right indicates that the system response is notably intense around some sections of the 
horizontal axis. This strong response arises when the kick frequency coincides or is close 
to one of the natural frequencies of the condensate, which are given by uj\ and UJ2 in 
equation (|6.35[) |94j . This effect is analogous to the resonances of mechanical systems, 
which occur when the driving frequency coincides with one of the system's natural 
frequencies. In this way, the strong resonance driving triggers particle dissipation and 
hence destroys the condensate. We point out that this is different from instability due to 
linear resonances, which take place when the contributions from different kicks interact 
coherently. In fact, these latter resonances are known to lead to chaos in the classical 
realm [S5J H3 M E3 ■ 

Here we have studied the wave function of a double-kicked BEC as well as the stability 
of the system. We have derived analytical expressions for the condensate wave function 
in the linear response regime. From the solution of the Gross-Pitaevskii equation, we 
identified resonance conditions for which particle dissipation is maximum. Our findings 
were verified by numerical simulations using a fast Fourier transform algorithm. We also 
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Figure 6.3: Left. /3 = 1.96, e = 0.04. In the strong repulsion regime NNP grows 
exponentially and the condensate becomes unstable. Time is measured in numbers of 
pairs of kicks 

. Right. Slope b characterizes the intensity of particle leakage in the condensate. As 
can be seen, diffusion peaks around some values of g, evidencing the effect of 

resonances. 



investigated the second quantized Hamiltonian and established useful formulae regarding 
the dynamics of the NNP. Our study is relevant to experiments of cold atoms in optical 
lattices and our results are interesting from a purely theoretical viewpoint. 



6.4 Advancements in the field 

Since the presentation of our results in 2007, there has been further progress on the 
subject of condensate stability. Here we want to summarize the main results of these 
interesting works in order to provide a complete view of the actual state of the topic. 
A complete characterization of the condensate stability has been reported in [84] . in- 
cluding a map of the stability parameters showing the position of the resonances for 
weak and strong interactions. In such reference the authors study how the non-linear 
response due to the interaction among atoms affects the well understood linear response, 
for which linear resonances have been studied and characterized. In this latter case, it 
is known that resonances are located at integer multiples of the so called Talbot time, 
T = Air. It was found that the non-linear dynamics not only displaces the positions 
of the linear resonances but also affects the nature of the response. Specifically, it is 
shown that in the range 1 < g < 20 the resonance response is characterized by a sharp 
cut-off of the NNP and therefore the condensate revives after it has lost stability. The 
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non- linear resonances manifest individually for small values of g and k, but they pro- 
gressively proliferate and overlap as the parameters are slowly tuned up. Likewise, the 
results of simulations using the GP equation are compared against the results of a sec- 
ond quantized model constructed from the insight provided by previous studies of the 
kicked-rotor dynamics [86j . It is shown that by including Bcliaev and Landau terms 
coming from phonon-phonon interactions it is possible to reproduce the sharp fall in 
NNP close to a resonance as well as the shift in the position of the resonance in regimes 
where depletion of the ground state is < 10%. Conversely, the same approach fails to 
reproduce the condensate dynamics when such terms are not taken into consideration. 
The second quantized model clarifies many issues regarding the depletion process and 
helps understand how resonances overlap and also how the cut-off process takes place. 
In the second quantized approach the condensate revival appears as a combination of 
Beliaev and Landau processes, in a way that is very similar to non-linear self-trapping 
of a BEC in a double-well potential. Certainly, the cut-off phenomenon can be under- 
stood by the following first principle analysis. When the condensate is being kicked 
slightly off-resonance, the interaction drives particles from the condensate towards the 
first excited mode. This triggers particle migration to the second mode, which at the 
same time induces a phase shift that brings the first mode closer into resonance. This 
is known as non-linear feedback and it is specially responsible for the strong condensate 
depletion seen before the cut-off. However, if the kicking is not sufficiently close to res- 
onance, the synchronization process among excited modes does not occur and depletion 
is suppressed, which results in the abrupt cut-off reported in reference [84] . Further 
inspection suggests that the intensity of the cut-off is much sharper than in double- well 
condensates. Moreover, the resonances are classified in two basic categories. First, linear 
resonances, which derive from the Talbot resonances and are seen to be paramctrically 
displaced by the non-linearity. These lead to strong, high amplitude response. Second, 
non-linear resonances, which disappear when g = 0. Such resonances can arise individu- 
ally, as Bogoliubov modes, or in combination of two or more modes. Non-linear response 
is weaker than linear response, but surprisingly, a Lyapunov analysis reveals that non- 
linear resonances are unstable, contrary to linear resonances. A possible explanation of 
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this is that non-linear resonances are linked to exponential oscillations and therefore the 
Lyapunov stability criterion becomes questionable. Finally, the stability profile could be 
used in applications where very sharp excitation thresholds are needed to carry out high 
accuracy measurements, as for example measuring gravity or detecting small changes of 
frequency in rotating BEC's. 
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Conclusions 



In this work we have adjusted recently developed numerical methods to probe the ground 
state as well as the dynamics of bosonic systems. In addition, we have introduced a 
numerical method that displays some advantages over standard approaches, and we 
applied it to solve challenging problems. 

In the first part of the document, we presented the concepts and ideas that we subse- 
quently used along most of the work. As numerics is an important part of this research, 
we extensively discussed how our programs were constructed and how they work. We 
then discussed alternative uses of MPS in cases where analytical solutions for the Hcisen- 
berg operators are available. The efficiency and reliability of our method was tested, 
first in simple scenarios, and then in more complex situations. The most important ap- 
plication of the method in this work was to find the state of a system of bosons initially 
prepared in a separable state after a many particle collision has taken place in the centre 
of the chain. This calculation would have been much more difficult to do using TEBD or 
tDMRG as both of these methods are suitable only for short-time intervals. Likewise, we 
also derived analytically the conditions for an efficient collection of bosons on the chain 
terminals after the particle waves collide. This is an important point as in the absence 
of controlling mechanisms boson waves flatten and entanglement diffuses. Conversely, in 
our approach particles are collected so that the entanglement contained in the system as 
a result of the collision can be used. Our graphs show that the entanglement generated 
in this way is quite substantial and it is shown that such entanglement can be made to 
grow by just adding more particles to the process. 

We next focus on applications of TEBD to the BH model, aiming at enhancing the 
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amount of entanglement between the ends of the chain, either in the ground state or 
as a result of dynamics. We show that EEE in the ground state scales logarithmically 
with a positive coefficient when a hopping profile with perfect transmission properties 
is incorporated in the Hamiltonian. We conclude that such effect is due to a marked 
change in the tunnelling profile, which undergoes a transition from local to global scope. 
We use the real-time version of TEBD to simulate the dynamics of a chain with high 
repulsion in intermediate sites. The amount of entanglement generated out of dynamics 
alone is smaller than what can be obtained from the ground state. We then use a 
scheme that combines both ground state and dynamics to generate entanglement. Our 
results indicate that it is in fact possible to use dynamics to increase the amount of 
entanglement contained in the ground state by taking advantage a pcrturbative scheme 
that we have introduced. 

In the last part we concentrate on kicked bosonic systems. We utilized a perturbative 
approach to obtain an analytical expression for the bosonic cloud in the linear response 
regime. Additionally, wc derived useful formulae regarding the number of non-condensed 
particles. We find that the condensate is highly depleted by a kicking with a driving 
frequency that matches any of the natural frequencies of the system. Similarly, we have 
reviewed the latest advances in the field. 

Among the several potential research extensions of the present work we would like to 
mention the possibility of studying chains with periodic boundary conditions, where 
vortices are expected to appear. Also, the numerical method presented here is liable of 
improvement. Indeed, it seems plausible that the method can be formulated without 
the need for an explicit solution of the Heisenberg equations of motion for the opera- 
tors. The procedure can be applied to a wide variety of problems, not only in linear 
arrangements, but also in lattices, since in our method the geometry of the system does 
not interfere with the reduction process. Potential extensions of the method to spin 
and fermionic systems are quite direct. Additionally, it would be very interesting to 
implement the configuration proposed in chapter [5] in an actual experiment, since such 
configuration can be realized in optical lattices using state-of-the-art technology. One 
of the main conclusions of this work is the enhancement of entanglement in chains with 
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PTH hopping. Similarly, as repulsion constants are turned up in the middle, we can 
think of intermediate places as forming a quantum wire. Therefore, it should be inter- 
esting to study a system made up of two particle-reservoirs connected by a quantum 
wire. As the infinity-repulsion BH model is known to derive in the XX model, it is 
likely that a problem like this can be solved analytically. If this is the case, it should be 
possible to establish the actual relation between the hopping profile and entanglement. 
The present has been an interdisciplinary investigation that covered areas such as numer- 
ical methods, many-body systems, quantum information and so on. Our findings help 
understand how entanglement relates to some physical processes that occur in bosonic 
systems. Such understanding can be used in different ways, as for example, to propiti- 
ate the circumstances that are convenient to produce entanglement in an experiment, 
or just as a way of appreciating the phenomenology of bosonic systems from a different 
perspective. Altogether, we hope our contribution to be sufficiently relevant so that it 
can inspire further investigation in any of the very exciting scientific fields with which 
this research overlaps. 
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Appendix A: Publications 



• Jose Reslen and Sougato Bosc, Long Range Free Bosonic Models in Block 
Decimation Notation: Applications and Entanglement, arXiv:0907.4315 
Submitted. 

• Jose Reslen and Sougato Bosc, End-to-end entanglement in Bose-Hubbard 
chains. Physical Review A, 80: 012330, (2009). 

• J. Reslen, C.E. Creffield and T.S. Monteiro, Dynamical instability in kicked 
Bose-Einstein condensates: Bogoliubov resonances. Physical Review A, 

77: 043621, (2008). 
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